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A COMPUTATIONAL METHOD FOR TIME-OPTIMAL SPACE RENDEZVOUS 


By Ernest S. Armstrong and Athena T. Markos 
Langley Research Center 

SUMMARY 

A method is presented for obtaining space- rendezvous trajectories between a 
one-stage rocket vehicle and a target in a circular Keplerian orbit. The trajectories 
satisfy the necessary conditions of the Pontryagin maximum principle for time- optimal 
rendezvous in which no terminal mass constraint is placed on the rocket. The use of 
Pontryagin T s theory leads to a two-point boundary- value problem. A digital program is 
given for the iterative solution of this problem. The method is successfully applied for 
the determination of time-optimal rendezvous trajectories between a vehicle launched 
from the surface of the moon and a target in an 80-nautical-mile circular orbit. 

INTRODUCTION 

The study of trajectories for rendezvous in space has been of interest for many 
years. An early summary of some aspects of the problem is given in reference 1. The 
application of optimization theory to trajectory computation has also received attention 
(refs. 2 and 3). Paiewonsky and Woodrow (ref. 4) have considered the problem of a 
three-dimensional time-optimal space rendezvous between a single-stage rocket vehicle, 
with constrained terminal mass, and a target in a circular Keplerian orbit. Linearized 
dynamic equations were used to describe the motion of the vehicle, and the mass limita- 
tion was of the form of a terminal inequality constraint. Linearization of the equations 
imposes the condition that the rocket be in proximity to the target. Reference 5 has 
considered the three-dimensional time-optimal rendezvous problem with unsimplified 
dynamic equations and fixed terminal mass through the dual problem of three-dimensional 
fuel-optimal rendezvous with specified final time. The present paper continues the 
extention of the problem of reference 4 to unsimplified dynamics by considering three- 
dimensional time-optimal rendezvous with unspecified terminal mass. 

The mathematical model of reference 5, which treats the rocket as a point mass 
and takes into account rotation of the attracting center, is employed. The Pontryagin 
maximum principle (ref. 6) is applied to find the correct thrust magnitude and direction 
for time optimality. This operation leads to a two-point boundary-value problem in which 
certain initial conditions on a set of differential equations introduced by the maximum 


principle have to be found such that certain terminal conditions are met. Following a 
method developed in reference 5, a digital program is written to solve the boundary-value 
problem by iteration. The procedure is illustrated numerically by solving the problem 
of finding time-optimal trajectories of a rocket vehicle launched from the moon to rendez- 
vous with a target in an 80 -nautical- mile circular orbit. 

In addition to extending the work of reference 4, the digital program and accom- 
panying analysis provide a useful method for obtaining time-optimal rendezvous trajec- 
tories and control laws. Therefore, the digital program is discussed and a listing 
included. 


SYMBOLS 


A,M,N,K,L,Tj,T 2 constant matrices 


A i = /3 




B 


B ik = 


X 7 # 

seven-dimensional diagonal matrix with elements bi 


*7 (#) 3 


bj positive weighting elements (i = 1, 2, ... 7) 

3n 2 R s 3 Vi v k 


c ik = 


H 5 

effective exhaust velocity 


Dik 


3f2 2 R s 3 (2i^ i v k + d) 

(*) 5 


d = !// 2 ( Xl + Rs x ) + ^ 4 (x 3 + R Sy ) + ^g(x 5 + Rs z ) 


Eik = 


15S2 2 R s 3 dViV k 


(k)’ 


EleGi)] scalar measure of terminal error, 


’ btfe^Sj] 2 


l 

i=l 


2 



e 


seven- dimensional vector with elements e^fa) 


e^ error criteria (i = 1, 2, ... 7) 

ik _ 3 ^ 2R s 3 (^ e v i + ^m v k) 

em " " (*f 


H 

I 

i,j,n 

i,j,k 

M 

M 

m 


R s 

R s 

R S X ? R Sy> 


pseudo-Hamiltonian function of the maximum principle 
identity matrix 

integers (i = 0, 1, . . . n; j = 1, 2, ... 6) 
unit vectors 

constant matrix defined in appendix A 
maximum value of H with respect to u 
mass of launch vehicle 
initial mass of launch vehicle 
magnitude of R s 

vector from center of attracting body to target 

R s elements of R s 

z *”* 


R s first derivative of R s 

R s x ’ R s y > R s z * irst derivatives of R Sx , R Sy , and R Sz 
Ry magnitude of R v 

R v vector from center of attracting body to vehicle 

Ry x ,Rv ,Ry z elements of Ry 


1 1 uni 


t 


Rx; 


first derivative of R v 


R v x ? R Vy^v z 
r = R v - Rg 


first derivatives of R Vx , R v , and R v 


r x? ry,r z elements of r 


r x ,ry,r z first derivatives of r x , ry, and r z 


dummy integration variable 


1 if p > 0 

sgn p signum function defined by (-1 if p < 0 

Unspecified if p = 0 

T magnitude of thrust vector 


thrust vector 


time 


initial time 


[*2^3] 

t' 


Ui 


final time 

nonzero subinterval of |”t O’tf] 
isolated point at which M’i//(t) = 0 
elements of u (i = 1, 2, 3) 


u 4 


magnitude of T 


u = u 4 u 


u 


unit vector in direction of u 


u. 


optimal form of u^ 


H im III 





1 


- * * -*■ * 
u =u 4 u 


optimal form of u 


v i’ v k 


variable v^, v 3> or vg 


V 1 = X 1 + R sx 


v 3 - Xg + R S y 


v 5 ~ x 5 + R s z 


x, y,z rotating axis system defined by figure A-l 


x ,y ,z 
X.Y.Z 


inertial axis system defined by figure A- 2 


^ = yfa + R Sx ) + (x 3 + R Sy ) 2 + (x 5 + R Sz ' 


X = COl^Xj, 


state variables (i = 0, 1, . . .8) 


first derivative of x. 


initial values of x^ 


first derivative of x 


initial value of x 


Y{x,x^j vector defined in equation (A9) 

ai,a j unknown parameters 


seven-dimensional vector with elements ai 


bound on thrust magnitude 


5 


5a 


5a 


7 

®c> 


M 


correction to a 


variation in a. 


control angles (see fig. A- 3) 


e v °,(p v ° initial vehicle angles (see fig. A-4) 
^o^o’^o angles determining target orbital plane 
X parameter governing step size of 5a 


universal gravitational constant multiplied by mass of attracting body 


switching function 




variables introduced by the maximum principle, with the sub- 
scripts equal to 0, 1, ... 8 


tf'i 


first derivative of 


if/ = coi^j, . . . i/'g) 

xf/ first derivative of i// 

+ ^4 2 + ^e 2 

Cl constant angular velocity of the target in its orbital plane 

( x> angular velocity of body about axis of rotation 

Mathematical notation (with arbitrary symbols used as examples): 

a(t) first derivative of a(t) with respect to t 

a(t) second derivative of a(t) with respect to t 

n 

a • b = ^ ajbj where a = col^aj, . . . a n ) and b = col^bj, . . . b n ) 
i=l 
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||al| = (a - a) 1/2 


[a,b] 

[a,b) 

9b (a) 
9a 

8 bj(a) 

9aj 


e 

i 


0 

-1 


closed interval 

interval closed at a and open 
Jacobian matrix with elements 


at b 

9bj(a) 
Cij= 9a j 


first partial derivative of bi(a) with respect to aj evaluated at 
a * col(ai, . . . a n ) 


belongs to a set 
denotes matrix transpose 
null vector 

as superscript to matrix, denotes inverse 


DEVELOPMENT OF THE BOUNDARY- VALUE PROBLEM 


In appendix A the dynamic equations are developed for a space vehicle which is 
required to rendezvous with a target in a circular Keplerian orbit about a rotating body. 
The vehicle is a one- stage rocket, treated as a point mass, with bounded thrust magnitude. 
The controls are the magnitude and direction of the thrust vector. 

From equation (A8), the dynamic equations when written as first-order differential 
equations in relative coordinates take the form 


*0 = 1 

£ 

11 

i u 4 Mu 

x = x 7 + Y ( x ’ x s) 

(x(t 0 ) = x 0 ; x(tf) = o) 

k --!H 

x 7 c 

J? 

o, 

II 

3 

o 

*8 = 1 

(x 8 (to) = to) J 
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The vector x is a six-dimensional vector whose elements Xj, Xg, . . . Xg are 
the relative position ^x 4 ,Xg,x^ and velocity (xgjX^Xg) components of the vehicle and tar- 
get. The variables x 7 and Xg are the instantaneous vehicle mass and the time, 
respectively. The vector Y(x,Xg), given by equation (A9), includes the kinematic parts 
of the equations. The scalar u 4 and unit vector u are, respectively, the magnitude 
and direction of the thrust vector T. The matrix M is a constant matrix defined in 
appendix A. 

Pontryagin's maximum principle (ref. 6) is now applied to find controls which mini- 

ptf 

mize \ dt while satisfying equation (1). In using the maximum principle, a new vari- 
to 

able Xq’ is introduced such that xq = l(xg(t 0 ) = o) and the problem becomes one of 
minimizing xq ( tf). From reference 6 the following conditions must be satisfied: 

(1) A function U4 ^ p and a function u ^with ||u|] = 1^ must be chosen to 
maximize 

H^ 0 , • • • ^ 8 ; x p . . . x g ; u p . . . u^ = ^ 0 x Q (t) + xp( t) • x(t) + i// ? (t)x 7 (t) + ^ 8 (t)x g (t) 

( 2 ) 

for fixed i//. and Xf (i = 0, 1, . . . 8). The vector \p is equal to col^j, • • • i/'gj. 
The terms (i = 0, 1, . . . 8) are nine additional variables introduced by the 

Pontryagin maximum principle and defined by 

*>. «) 


(2) For any te£t 0 ,tfj, i//g(t) = Constant = 0 and the maximum value of 

H(i// 0 , . . . ^ 8 ; Xj, . . . x 8 ; Uj, . . . u 4 ) with respect to Uj, . . . u 4 , given by 
M^i }/q, . . . i//gj Xj, . . . XgJ, must be identically zero over jt 0 ,tfj. 

(3) The transversality condition must be satisfied. 

Since x(tf) is specified and x 7 (tf) and Xg(tf) are unspecified, the transversality 
condition discussed in reference 6 yields 

^ 7 (ff) = i^g (tf) = 0 

From a consideration of condition (1), 


H = ^0 + ^ * 

whereby, from equation (3), 


u 4 Mu 


+ Y| 


: (i,x 8 ) 


V4 


+ 


8 





^0 = ° 


^ = 


^7 = 




3x 

U^l// * Mu 


X n 


®Y(x,Xo) _ 


(*0 = °) 

(^(to) undetermined) 

(W = °) 

(^ 8 (tf) = 0) 


> ( 4 ) 


If M> # 0, the u which maximizes H and satisfies ||u||=l is 

_ MV 
"||MV|| 

since i// ■ Mu can be written as MV * u. Then H becomes 

s = ^0 + " 1 ?) + ^ ' ? (*’ x 8 ) + ^8 

and the function u^ = /3 which maximizes H, if x 


(5) 


— does not vanish identi- 


cally over a nonzero interval in [t 0 ,tfj, is 

u 4* = |(! + s e n p) = v 


(p(t) > 0)} 
(p(t) < 0)/ 


( 6 ) 


with 


_ HmvII 'h 


The complete control now takes the form 


u* = u 4 *u* = |(1 + sgn p) M jl- 

2 ||MV|| 


(7) 


The function p(t) is the switching function for the system. The function p(t) has 
no zeros on (t 0 ,tfj except possibly at tf, which simply means the thrust is always at its 
maximum value. Since ^(t) is such that 


= |(1 + sgn p) 2 


Mj^/ ^ Q 


9 



/ 


■ I ■ I ■ I ■■ I •• 


and ^(tf) = 0, it follows that xf/^( t) S 0 for all te|t 0 ,tfj. If ^ 7 (t') =0 at some 

t’e[to, tf), the condition p(t) < 0 must be satisfied for t > t' in order to meet the con- 
dition ^(tf) = 0. The implication is that from t' until tf, the vehicle coasts; that is, 
T = 0 from t' to tf. However, to rendezvous at tf requires that the vehicle and 
target have the same position and velocity. Therefore, ^(t) can vanish only at tf, 
since it is not possible for the vehicle to coast into the same position and velocity as the 
target; that is, i//^(t) < 0 and p(t) > 0 for te[t 0 ,tf). Equation (7) then takes the form 


M'^l 


(8) 


and M becomes 


M= + /3p + ^ • Y(x,x g ) + 

In an optimal system, M'^ = (^2’^4’^6'j = ® over a nonzero interval (for example, 

over po>tf|^ cannot be allowed. For, if it is, the differential equations for 

\p 2 > and i//g imply = 0 over ^t^t,^. In fact, xp(t) = 0 would be the 

solution of the \p system over |t2,tfj. Since 1/^7 (tf) and iZ'g(tf) both vanish, ^(t) s 0 

and i//g(t) s 0 over * Also ¥ - 0 implies that = 0 over {t2,tfj, giving 

(* 0 ,*(t),* 7 (t),* 8 (t)) = 0 over [t2,tfj, which is a contradiction to the maximum principle. 

The maximum principle states that for each te [to, tf], the vector ^Q(t),i//j(t), . . . i//g(t)) 
is nonzero. Isolated points at which M'tp = 0 have no effect on the solution, since u(t) 
is bounded. For definiteness, 

i(f) = lim 

t— t - ||m’ 

at any isolated points t’ where M'\p(t') = 0. 

Equations (1) and (4) now take the form 

*0 = 1 


x = 


/3MM* i// ^ 2 R s 3 


M’ 


Ax + R 


s ( x 8)ir 


fopo) = 0) 

[Nx + MRs^Xgjj + n 2 MR s (x 8 ) 




+ (n* + 2wK + w^L^x 
= _£ 


k 8 = 1 


(x(t Q ) = x 0 ; x(tf) = 0 ) 
(x 7 (t 0 ) = m 0 ) 

(Xg(t 0 )= t Q j 




> (9 a) 


10 



and 




* 0 = ° 


(*0 - 0)1 


fi 2 R,, 3 N , ii/ 
* = S 


3J2 2 R s 3 


<3 c' (|nx + MR s (x ft )j • iAa’[Ax + RsfxjAj 

ax + r s (x 8 )|| 3 ||ax + r s (x 8 )|| 5 U- sU ^) J L sU ^1 


- (n + 2a)K' + w 2 L*)^ 


; gj[Mwj[ 

x,2 

S1 2 R • Mjt.fxo} 
] Ax + RftfXg 


^i//(t 0 ) undetermined) 
(^ 7 (tf)=0) 


> Ob) 


A,-* ||5 '( R s( x 8) * + R s( x 8)]} • [Nx + MR s (Xg)jj (V 8 (tf) = 0) 

J 


By comparing equations (9a) and (9b), a two-point boundary -value problem is recognized; 
\pQ, i//(t 0 ), ^ry(to), >^g(to)> and tf need to be found such that x, M, and i^g are 
zero at tf. This boundary-value problem can be placed in the simpler form 

M(^ 0 , • • • P Q ; Xj, x 8 ) = 0 

which yields 

* 8 (to) = -[* 0 + ^ (to) + Wo) • Y(x 0 ,t 0 )] (10) 

Since M = 0 over (t 0 ,tfj, since ^(tf) = 0, since x(tf) =0, and since Y^0,tf) = 0, 
it follows that 

M= p Q + J3p(tf) = 0 
or 

* 0 = -0P(tf) (ID 

The boundary-value problem then reduces to satisfying equations (9a), (9b), (10), and (11) 
and finding parameters i//j(t Q ), . . . i//^(t 0 ) and tf such that x(tf) and (tf ) are 

zero. There appear to be eight parameters and seven boundary conditions. However, 
because of the homogeneous form of their differential equations, one of the parameters 
<^f(to) (i = 1, 2, . . .7) can be removed by normalization — that is, by fixing its value. 
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V 

Also, by observing the differential equation for xp^, it can be noted that the parameter 
t 0 ) could be determined after the other parameters are found by setting 



The minimal combination is then six parameters (^(t 0 ) and tf, with one of the elements 
of p(to) normalized) and six boundary conditions (x(tf) = 0) . However, since the correct 
algebraic sign of an element of ^(t 0 ) may not be known a priori, the problem is best 
solved by finding seven parameters (xp(t 0 ) and tf, with (t Q ) normalized) such that 
the seven boundary conditions ^x(tf) = 0 and iprj (tf ) = 0) are met. 

SOLUTION OF THE BOUNDARY- VALUE PROBLEM 

Development of Iterative Logic 

The approach taken to obtain a solution of the foregoing boundary -value problem is 
that discussed in reference 5. In reference 5, given a similar boundary -value problem, 
a vector e(a) is defined such that when e(cF) = 0 the boundary conditions are satisfied 
and the unknown parameters for the problem are a. 

The vectors a and e(a) correspond to col(Q!x) and col(ei(a)) 

(i = 1, 2, . . . 7), respectively, with elements 

«! = *i(to) 
a 2 = 

“6 = 

«7 = tf 

The quantities Xf (i = 1, 2, . . . 6) and iprj are written Xf(o) and xf/, j(a) to indicate 
their implicit dependence on a. 

The magnitude of e(a) is measured by a scalar quantity 

2 

where B is a positive definite diagonal matrix of weighting elements. Here 

E [e(30] = |[b 1 x 1 2(S) + b 2 x 2 2 ( 5 ) + . . . + b 6 Xg2 (a) + b^ 7 2 (B)] (13) 

where bf > 0 (i = 1, 2, ... 7). 
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A value of a is assumed, the differential equations are integrated forward in 
time until t = tf = a 7 , and E[e(Z?)| is evaluated. If e(a) vanishes or, for practical 
purposes, is sufficiently small, the boundary-value problem is considered solved. Other- 
wise, the assumed a is corrected by 


5a = - 


B + XlT B e(a) 

_ da da J 9 « 


(14) 


Where X > 0 is adjusted such that 

Eji(cr + 6o)j < Eje(o)j 


(15) 


For this problem 


is given as the partitioned matrix 
9a 



! 



i 

9x(a) ] 

x(a) 

9e(<7) _ 



9a 

i 

i 

d^ia) 

da 

^ 7 (a) 


, 



(16) 


where 


is a Jacobian matrix with elements 
9a 


9i f/rj (a) 

j = 1, 2, ... 6) and = — is a row vector with 


9xi(o) 


9a 


elements 


(i = 1, 2, 

9j// 7 (o) 


9ai 


The first six elements of 6a" are the corrections on i^(to) 
The seventh element 6 a^ is the correction on the last value of tf. 
is given by a 7 + 6a 7> 


■ (j = 1, 2, . . . 6). 

(i = 1, 2, . . . 6). 
The next final time 


The procedure is designed to be applied iteratively and generate a monotone 
decreasing sequence of E[e(a)] converging to the smallest Eji(o)] available relative 
to the initial choice of a. Success of the method is dependent on the user's ability to 
find starting values of a and at each stage find values of X such that equation (15) is 
satisfied. At each stage a one-dimensional search may be performed to find the values 
of X. However, for this problem it was found that a value of X of 10 or 1 would suffice 
throughout. In general, the method does not guarantee a solution for an arbitrary 
boundary-value problem. It has, however, been highly successful in yielding solutions 
to the boundary -value problem under consideration herein and to others (see ref. 5). 
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Expanded versions of equations (9a) and (9b), with Xg replaced by t, are 


*i = x 1+1 (i = X, 3, 

5; Xi(t 0 ) = x io ; Xi (tf) = 0) 

» 2 2 3 *1 +Rs x 2 2 

x 9 - — _ ^Rg* 5 „ x + fTR Sv + a/xi + 20JX. 

2 x 7 # (vK) s X 1 4 

( x 2 (t 0 ) = x 2q ; x 2 (tf)=0) 

. 2 3 X 3 + *^ s v 0 9 

x 4 = — : - 0 R s — *-*- + J2 Ro - 2u>x 0 + 

4 X 7 # (vS) 3 Sy 2 3 

( x 4 (to) = x 4 0 ; x 4 (tf) = 0) 

' 

2 3 X 5 + R s z 2 

x 6 = — ^ - fTRg ^ + 0 2 R S 

6 x 7 # (^c) 3 

( x 6 (t 0 )=X6 o ; x 6 (tf)=o) 

x 7 = '§ 

(x 7 (t 0 ) = m 0 ^ 

and 



. fi 2 R s fy 2 2 3 Xl + R s x o 

^2 = '^1 + 2aM ^4 


^ 2 R S 3 ^4 2 3 X3 + Rs V 2 


(0 


^4 = - 2 o >^ 2 - 

ft 2 R s 3 iK „ ^ x 5 + r s 

_ S 0 on^D 


= —nssr - 3n R s d 


(v5F) 


(&Y 


h = "^5 



^l//^(to) - 0! 

(^ 2 (to) = a 2 ) 


(^3(to)=a 3 ) 
(^4 (to) =Qf 4 ) > 


( 17 b) 


(^ 5 (to)=a 5 ) 

(^ 6 (to)=« 6 ) 

(^7 (to) normalized; 1^7 (tf ) = 0) 

J 
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where 


and 


# = \/^2 2 + + * + ^6 2 


^ = ^/( x l + Rs xf + ( x 3 + R sy) 2 + ( x 5 + R s Z ) 5 


d = ^ 2 (Xj + R Sx ) + ^ 4 (x 3 + R Sy ) + ^ 6 (x 5 + R Sz ) 

The derivatives needed to form equation (15) can be obtained (ref. 7) by solving the fol- 
lowing system in conjunction with equations (17a) and (17b) from t Q to tf, with 
j = 1, 2, ... 6: 


d 



9x i+l 


dt\ dot 


d r x 2\ 


da 


J 

a^ 2 


/^(to) \ 

[-kr B °i i = 1 > 3 ’ 5 




at l toTj = A2 + ° 24 


+ Be 


3 


9Xr 


8 *4 9^ 

+ B 2 6 “ — + 


9a 


] 


9x 

+ C I o - + 2 co ** + C 1 c — — 

9aj 9aj 15 9aj 


9a j 


9x c 


'll 


q 2 -d 3 
11 Rs + „2 


(0 


9Xj 
9a -i 


d i Bx 4\ 8 ^2 8 ^4 8x 1 8x o 

-7- — - = B04 — — - + A 4 — - + Bar ^ + C n — i - 2co 


V 


9a i 


C 3 3- 


S2 2 R s 3 «, 

— + CO 2 


9a \ 


(«) s 


9x« 9 x r 

4 +C35 4 


d r x 6\ 9 ^2 8 ^4 9 ^6 9x l 9x 3 

dtUaj j = B 26 8^ + B 46 8^7 + A 6 8^ + c 15 g^T + c 35 g£7 


c 55 


S7 2 R s 3 

9x 5 

(V£) 3 

9aj 


^ x 2 (t 0 ) 


9a -i 


= 0 


> (18a) 


/, 9x 4 (t 0 ) 


9a i 


= 0 


^9Xg(t 0 ) 


9ai 


= 0 


J 


15 


y 


d/^A 


dtUo 


h 


Sl 2 R s 3 

w 


3^2 ^4 Sfe 8xi 

aofj 13 15 9q ,_. ( 21 ll)ao/j 




^(toj ri o = i^ 
Lo (i * i)y 


90H 


d ( d J% _ 2w _^4 _ 8 ^1 

dty30!jy 8(3! j 8(3! j 


_ r 1 (j=2 ^ 


3a n - 


J> (j * 2 )J 


i^ = -c 13 ^ + 


dtiaa 


3qj 


3 


« 2 «S 3 c 2 

- C 33 - CO 


w 

8x 0 


9^4 


So;.; 


(°43 + E 33)^ + ( F 46 + E 35)a^ 
9^2 d ^3 


Cq *1^ + (f 13 + E ^ 
8o!i V 42 13/ dot i 



±L = -2co 


^^ 3 (t 0 ) _ ri 0 = 3)\ 

LO (j * 3 )/ 


8(3! i 


dt\8a 


8(3! j 9(3! j 


'8^4 (to) _ ri (j = 4)> 

Lo (j * 4 )y 


da 


d /8^\ 8^ 2 9^4 

dtfer) = ' Cl5 — - c 35^ + 


8or i 


8(3! i 


s/ 2 r s 3 




c 55 


91//, 


J 

Bx* 


/ O C \ BXp- 


f + (f|| + e 15 )^ 

/ 9^ 5 (t 0 ) ri a = 5 f\ 


da i 


J> (j#5)/ 



dt \ 80 ! 


9l/ / c 

t 

Sqm 


' 9 ^ 6 (t c ) _ f 1 (3 = O) 1 ) 
Lo (j * 6V 


8(3!; 


dYf^ 

dt v“i/ 


/ 8^ 8 xI/a 9i//g 

2 8^ + ^4 9^ + ^6 8^ 


X, 


2^ 


In this system 


Ai = /3 




|x 7 # x 7 (#) 3 J 


% 7 (t 0 ) 


Ba-i 


= 0 




(18b) 
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and 


with 


and 


B ik = 


C ik = 


_ j^k 
x 7 (#) 3 

3Q^R s ^VjV^ 

w 



15J2^R s ^dvjVj c 

(n/5c) 7 


P ik 

em 


3C2 Z R 


2-r 3 


(*■ 


e v i + V v k 


(^) E 


V 1 = X 1 + R s x 

v 3 = x 3 + R sy 
v 5 = x 5 + R s z 


Iteration Sequence 

In summary, the procedure used to solve the boundary-value problem presented in 
the preceding section is as follows: 

(1) Assume a value of a given by equation (12). 

(2) Solve equations (17a) and (17b) to tf = a 7 and evaluate e(a) given by equa- 
tion ( 12 ). 

(3) If e(a) is sufficiently small, then the problem is solved; a gives the unknown 
initial conditions and final time, and the corresponding solution of equation (13) gives the 
optimal trajectory. 

(4) Otherwise, with the use of the solution given by equation (13), for the assumed 
a, solve equations (18a) and (18b) and evaluate equation (16). 

(5) If previous value is not satisfactory, find A such that equation (15) follows. 

( 6 ) Replace oi by a + 5a and return to step (2). 

A digital computer program for this procedure is discussed in the next section. 
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DESCRIPTION OF PROGRAM 
General 

The program was written in FORTRAN IV language for the Control Data series 6000 
computer systems at the Langley Research Center. A complete listing of the program is 
found in appendix B. 

Upon acceptance of an assumed value of a and appropriate system constants 
characterizing the particular rendezvous problem to be solved, the program proceeds 
according to the steps listed in the preceding section. The program does not contain a 
search routine for the determination of X in step (5). It is expected that in each 
instance a value of X which will work for the complete iteration process can be found. 

The mathematical symbols used in the theory and their FORTRAN equivalents are 
given in table 1. 

TABLE 1.- FORTRAN EQUIVALENTS OF MATHEMATICAL SYMBOLS 


| Mathematical symbol 



FORTRAN equivalent 

bi 

(i = 1, 2, . 

• ■ ?) 

B(I) 

(I - 1, 2, ... 7) 

c 



C 


ei(a) 

(i = 1, 2, . 

- • 7) 

E(J) 

(J = 1, 2, ... 7) 

E[e(a)] 



EDP 





10 


Rs 



RS 


Rs x » R s y ,Rs z 



RSV(I) 

(I = 1, 2, 3) 

Rs x ,Rs y ,Rs z 



DRSV(I) 

(I = 1, 2, 3) 

R v 



RV 





RVX, RVY, RV Z 


Rv x , R v y ,Rv z 



DRVX, DRVY, DRVZ 

t o 



VAR(l) 


tf 



TF 


Tl 



XT1(I,J) 

(1= 1, 2, 3; J = 1, 2, 3) 

T 2 



XT2(I,J) 

(1= h 2, 3; J = 1, 2, 3) 

u l» u 2> u 3 



UX,UY,UZ 


x i 

(i = 1. 2, . 

- 7) 

VAR(I) 

(I = 2, 3, ... 8) 

€-b4 + *i 

da 3a 



PEMAT(I,J) 

(1= 1, 2, ... 7; J = 1, 2, ... 7) 

0=1,2,.. 

3«j 

•6; j = 1, 2, . . 

. . 6) 

PX(I,J) 

(I * 1, 2, ... 6; J = 1, 2, ... 6) 

“i 

(i = 1, 2, . . 

• 7) 

PEVEC(1, 1) 

(I = 1, 2, ... 7) 

0 



BETA 


Oo 



THETAO 


0 V ° 



THETVO 


X 



LAMBDA 


M 



MU 





PHIO 


0 O 



PHrvo 


'Pi 

(i = 1, 2, . . 

. 7) 

VAR(I) 

(1= 9, 10. . . . 15) 

0) 



S0MEG 


a 



C0MEG 


8Cj 

.6; j = 1, 2, . . 

• 6) 

PPSI(I,J) 

(I = 1, 2, ... 6; J = 1, 2, ... 6) 

9a i 

(i = 1, 2, . . 

- 6) 

ppsn(i) 

(I = 1, 2, ... 6) 
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Subroutines 


Seven subroutines are used in addition to the main program. The purpose of each 
is given in table 2: 


TABLE 2.- SUBROUTINE LISTING 


Subroutine 

Purpose 

DERSUB 

Evaluates all differential equations to be solved 

CHSUB 

Tests for the end of a trajectory 

COMP 

Evaluates the position and rate of R s (t), the vector from the 


origin to the target 

ITERAT 

Computes and applies the correction 6a to the initial a 


ax. (5) aiMS) Ma) 

BLOCK DATA . . . 

Initializes — , — , and 

8 “j 8 «j 8a j 


(i = 1, 2, ... 6; j = 1, 2, ... 6) 

INT2 

Numerically integrates the differential equations with a 


fixed- step size method by employing a fourth- order 
Adams -Bashforth predictor formula and a fourth- order 
Adams- Moulton corrector formula 

MATINV 

Obtains the inverse of the matrix B + XI 

Ida da J 


Input 

Input is of the form shown in table 3: 

TABLE 3.- INPUT DATA 


Card 

FORTRAN 

Description 

FORTRAN 

number 

variable name 

format 

1 

NO 

Case number 

120 

2 

S0MEG, BETA, C, TF 

0), 0, C, tf = Ctrj 

4E20.8 

3 

PHIVO, THETVO, RV 

^V 0 ’ ^v°> Rv(to) 

3E20.8 

4 

DRVXO, DRVYO, DRVZO 

Rv x fto) » R-Vy (^o) 9 ^Vz (to) 

3E20.8 

5 

PHIO, THETAO, IO, RS 

o * @ 0 > L Of Rs 

4E20.8 

6 

VAR(l), VAR(8), MU 

to, -^0, M 

3E20.8 

7 

VAR(9) to VAR(12) 

a 1 to a 4 

4E20.8 

8 

VAR(13) to VAR(15) 

« 5 , a 6 , (to) 

3E20.8 

9 

Cl, SPEC 

Computing interval, printing fre- 
quency (see "Output" section) 

2E20.8 

10 

IPRINT, IER0R, IMAT 

(See "Output" section) 

3120 

11 

LAMBDA, CRIT, MAXIT 

X, stopping criterion, maximum 
number of iterations (see 
"Output" section) 

2E20.8, 120 

12 

B(l) to B{4) 

bi to b 4 

4E20.8 

13 

B(5) to B(7) 

bs to b 7 

3E20.8 
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All the input variables are dimensionalized and angles are in radians; a' i 

(i = 1, 2, . . . 6), bf (i = 1, 2, . . . 7), and (i = 0, 1, . . . 8) are considered 

dimensionless. 


Output 


The program offers several options for output. Regardless of the options, the input 
data are always printed initially. Afterwards, output is presented at each iteration 
according to the following input variables: SPEC, IPRINT, IER0R, IMAT. 


SPEC specifies how often results are to be printed. If SPEC = 10 1®, output will be 
printed only at t = to and t = tf. If SPEC = nCI, where Cl is the integration computing 
interval and n is a positive integer, results will be printed every n integration step. 

At t = t 0 , the variables that are printed are t 0 , i/q(to) (i = 1, 2, . . . 7), and uf 
(i = 1, 2, 3). At any other time t, determined by SPEC, the variables that are printed 
are t, i//j(t) (i = 1 , 2, . . . 7), R s (t), tg(t), R v (t), R v (t), xi (i = 1 , 2, . . . 7), 

Uj (i = 1, 2, 3), ^(t), x 4 (t), x 6 (t), /3, [|R v .(t)||, and the relative distances and 

velocities + Xg^ + and -^}Jx + Xg^ + x^. At t = tf, E^e(«)| and 5a are 


also printed. 


and 


IPRINT provides the option for printing the partial derivatives 

9 ^7 


9l//f(t) 9xf(t) 


8of 


da 


da 


(i = 1, 2, 


6; j = 1, 2, 


6). If IPRINT = 0, the partial derivatives are 


not printed; if IPRINT = 1, the partial derivatives are printed. 


IER0R provides the option for printing the truncation errors for the differential 
equations. If IER0R = 0, truncation errors for the differential equations are not printed; 
if IER0R = 1, the truncation errors are printed. 


IMAT provides the option for printing the matrix 


9e’_(o0 B 

9a 


9 e(a) 

9a 


+ XI 


, its inverse, 


and the product of the two. If IMAT = 0, the matrices are not printed; if IMAT = 1, the 
matrices are printed. 

The program terminates when convergence is reached ^El[e(a)] S CRIt) or when the 
maximum number of iterations (MAXIT) has been reached. 


EXAMPLE COMPUTATIONS 


The use of the foregoing procedures is demonstrated by the problem of a space 
vehicle launched from the surface of the moon to rendezvous, in a minimum of flight time, 
with a target in a circular orbit. The vehicle has a bounded thrust magnitude which, along 
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with the thrust direction, acts as a control variable. There is no terminal mass con- 
straint. The system constants used in this study are given in table 4. 




TABLE 4.- SYSTEM CONSTANTS 


Initial time, t Q , sec 

Upper bound on thrust magnitude, j8, lbm (kg) 

Initial mass, m Q , slugs (kg) 

Effective exhaust velocity, c, ft/sec (m/sec) 

R v (t 0 ) set equal to radius of moon, ft (km) 

Rg set equal to radius of 80-nautical-mile satellite 

circular orbit, ft (km) 

Universal gravitational constant multiplied by the 

moon mass, p, ft3/ sec ^ (m3/sec2) 

Angular velocity of moon about axis of rotation, u>, rad/sec 


0 

3504 (1589.4) 

285.5 (4166.3) 

9853.2 (279.86) 

. . . 5.707 X 10 6 (1739.4) 

. . 6.1934 X 106 (1887.7) 

1.727 X 10 14 (48.90 X 1011) 
2.66 x 10-6 


The satellite orbital plane was placed in the xy-plane of the rotating system 
(fig. A-l). Studies were made with the vehicle launched from rest from the surface of 
the moon with the launch site both in and out of this plane (planar and nonplanar case, 
respectively) . 

It was found that a workable set of bj (i = 1, 2, . . . 7) and X for convergence 
was bj = b0 = b§ = b7 = 1, b2 = b4 = bg = 10, and X = 10 or X = 1. These values were 
used throughout. It was observed that increasing the value of X yielded a slower con- 
verging process, while a decrease was apt to produce divergence. It was also observed 
that by increasing a particular bi, greater influence could be applied to the correction of 
the error e^; that is, this error would be corrected more quickly than before, but at the 
expense of the other errors. 

It was found that in the planar case, with the vehicle launched such that the satellite 
lead angle cp Q - cp v ° was 9° (with cp Q = 89°), the set of values 

a 1 = -100R Sx (t o ) = 7.8565 
a 2 = Rs x (to) = 5.296 x 10 3 
0*3 - Rsy(to) = -4.5016 
CL 4 = R Sy (to) = -92.44 
«5 = Rs z (to) = 0 
a g = Rs z ^o) = 0 
aq = 500 seconds 
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(ipijito) having been normalized at -1.184 x 10 5 ) leads to convergence with the values of 
bi previously mentioned and X = 10 in 12 iterations. This procedure gave a solution 
to be used as a nominal, or guessed, solution for neighboring lead angles. Table 5 shows 
the planar results obtained. For each lead angle the iteration was stopped when 

x(tf) • x(tf) + 2 (tf) S 1 

with \prj(t 0 ) normalized at -1.184 x 10 Since the results are planar, ^/ 5 (t) s j//g(t) = 0. 

It can be seen from table 5 that near the lead angle 13.7°, the smallest value of tf 
and hence the largest terminal mass occur. Figure 1 shows a graph of these results. 

TABLE 5.- TIME-OPTIMAL PLANAR RESULTS 


Lead angle, 

<Po ~ <Pv°, 
deg 

Unknown parameters 

<h(t 0 ) 


^3 (to) 

^4 (to) 

tf, sec 

8 

11.662 

5577.3 

5.6883 

2051.8 

547.9 

9 

12.929 

5948.3 

7.4703 

2638.8 

524.8 

10 

14.120 

6239.4 

10.026 

3422.5 

499.5 

12 

12.530 

5464.5 

17.992 

5427.7 

453.0 

13 

6.9429 

3774.3 

20.577 

5758.3 

443.3 

13.7 

2.4694 

2500.4 

20.507 

5501.0 

442.3 

14 

.73674 

2012.3 

20.108 

5315.8 

443.0 

16 

-6.2567 

-42.810 

15.897 

4005.4 

454.8 

18 

-8.4150 

-882.65 

12.256 

3110.3 

471.7 

22 

-8.5773 

-1472.0 

7.8560 

2155.4 

507.3 


Percent initial 
mass at tf 

31.7 

34.6 

37.8 

43.6 

44.8 

44.9 

44.8 
43.3 
41.2 

36.8 



Figure 1.- Percentage of initial mass at rendezvous as a 
function of lead angle <p 0 - <p^°. Planar case; 0 V ° = 0°. 
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Schematic views of the flight paths for planar solutions with lead angles 8°, 13.7°, 
and 22° are shown in figure 2. From this figure an idea can be gained as to the different 
maneuvers required for different lead angles. Arrows placed along the trajectories indi- 
cate the true direction of the thrust vector at 50-second intervals. The spatial coordi- 
nates used in this plot are to different scales. The xy-plane is viewed as being inertial 
since the total rotation of the moon was less than 0.1° for the longest flight time. 

Examples of out-of-plane results were obtained by holding cp Q and cp v ° fixed at 
93.7° and 80°, respectively, and allowing nonzero values of 9 V °. For 0 V ° = 2°, the 
planar solution for cp 0 - cp v ° = 13.7° was used as a nominal. The results for this non- 
planar case exemplify a typical sequence of iterations, and this sequence is tabulated in 
table 6. For 0 V ° = 2°, bj = b 3 = bs = b 7 = 1, b 2 = b 4 = bg = 10, and X = 1. Table 7 
shows other out-of-plane results for the fixed lead angle of 13.7°. A graph of the percent- 
age of initial mass at rendezvous as a function of out-of-plane angle 0 V ° is shown in 
figure 3. 



Figure 2.- Comparison of planar time-optical trajectories. <f\,° = 80° : <p 0 = 88°, 93.7°, and 102°; 0 V ° = 0°. 

(0.3048 meter = 1 foot) 
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TABLE 6.- TYPICAL TIME-OPTIMAL ITERATION SEQUENCE 
p 0 = 93.7°, cp v ° = 80°; 0 V ° = 2° (except for nominal results^ 


Iteration 


^l^o) 

^2 fro) 

Nominal 3 - 

2.4694 

2500.4 

1 

2.6048 

2526.5 

2 

3.5498 

2718.7 

3 

3.6100 

2732.8 

4 

3.6263 

2736.9 

5 

3.6325 

2738.6 

b 6 

3.6351 

2739.2 


Unknown parameters 


^3 (to) 

H (to) 

^5 (to) 

20.507 

5501.0 

0 

20.415 

5498.3 

-4.6177 

19.445 

5369.0 

-4.8561 

19.461 

5371.7 

-4.8560 

19.469 

5374.8 

-4.8587 

19.474 

5376.2 

-4.8596 

19.476 

5376.8 

-4.8599 


P 6 (to) 

tf, sec 

E (1(5)] 

0 

442.3 

1.75 x 10 10 

-1340.9 

443.2 

DO 

CO 

O 

X 

h-i 

O 

00 

-1409.4 

448.5 

2.23 x 10 6 

-1412.6 

448.8 

9.29 x 10 

-1413.6 

448.8 

9.55 

-1413.9 

448.8 

1.61 

-1414.0 

448.8 

.271 


a 0 v ° = 0; Xj(tf) = -582.7, x 2 (tf) = -0.03386, x 3 (tf) = -3275.5, x 4 (tf) = 0.15997, 

Xg(tf) = 1.8728 x 10 5 , x 6 (tf) = -27.581, and ^ 7 (tf) = 9.9489. 

b Xl (t f ) = -0.0017, x 2 (tf) = 0.230, x 3 (tf) = -0.006, x 4 (t f ) = -0.026, 

Xg(tf) = 0.001, xg(tf) = 0.015, and ^(tf) = -0.042 for a computing time of 23 seconds. 


TABLE 7.- TIME-OPTIMAL OUT-OF-PLANE RESULTS 


j^ 0 - 93.7°, cp v ° = 80°j 


6 V °, deg 



Unknown parameters 



Percent initial 


^2^o) 

^ 3 (t 0 ) 

^4 (to) 

^(to) 

^6 fro) 

tf, sec 

mass at tf 

0 

2.4694 

2500.4 

20.507 

5501.0 

0 

0 

442.3 

44.9 

2 

3.6351 

2739.2 

19.476 

5376.8 

-4.8599 

-1414.0 

448.8 

44.1 

4 

5.6559 

3154.0 

16.644 

4913.9 

-7.8427 

-2482.2 

466.8 

41.8 

6 

6.6424 

3323.6 

13.215 

4201.8 

-8.6545 

-3010.3 

492.2 

38.7 

8 

6.5482 

3215.7 

10.282 

3495.9 

-8.2675 

-3132.6 

520.0 

35.2 

10 

5.9679 

2970.9 

8.0805 

2914.0 

-7.4739 

-3048.8 

546.6 

31.9 
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Figure 3.- Percentage of initial mass at rendezvous as 
a function of out-of-plane angle S v ° for fixed lead 
angle <p 0 - of 13.7°. 


The computational time for obtaining both planar and nonplanar trajectories was 
less than 1 minute. 


CONCLUDING REMARKS 

A technique has been presented for obtaining three-dimensional rendezvous tra- 
jectories between a one- stage rocket vehicle and a target in a circular Keplerian orbit. 
The trajectories satisfy the necessary conditions of the Pontryagin maximum principle 
for time- optimal rendezvous in which no terminal mass constraint is placed on the rocket. 

The use of Pontryagin's theory leads to a two-point boundary -value problem. Cer- 
tain initial conditions on a set of differential equations introduced by the maximum prin- 
ciple had to be found such that certain boundary conditions were met. A digital program 
was given for the solution of this problem based on an iteration method. Given assumed 
values of the initial conditions, which do not yield rendezvous, the program attempts to 
correct these values in such a way that rendezvous is more closely attained. Iterative 
use of this procedure gives a sequence of trajectories converging to one yielding 
rendezvous. 


A. 
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The program was successfully applied to a problem in which a space vehicle was 
launched from the surface of the moon and required to rendezvous with a target in an 
80-nautical-mile circular orbit. Both planar and nonplanar trajectories were obtained 
with equal ease in less than 1 minute of computational time on the Control Data 
series 6000 computer systems. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., November 8, 1968, 

125-17-05-10-23. 
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FORMULATION OF DYNAMIC EQUATIONS 

The dynamic equations for a space vehicle which is required to rendezvous with a 
target in a circular Keplerian orbit about a rotating body are derived in reference 5. 

The formulation of these equations is summarized in this appendix. The vehicle is a 
one- stage rocket, treated as a point mass, with bounded thrust magnitude. The controls 
are the magnitude and direction of the thrust vector. 

Let x, y, and z be Cartesian coordinates of a rotating axis system located in 
the center of the body with the z-axis through the axis of rotation of the body. The geom- 
etry is represented in figure A-l. 

z 



The vector R s (t) is from the origin to the target, R v (t) is from the origin to the 
vehicle, and oo is the angular velocity of the body about its axis of rotation. The rela- 
tive distance between the target and vehicle is given by 

r(t) = R v (t) - R s (t) (Al) 

Since the target is assumed to be in a circular orbit, it moves in its orbital plane 
at a constant distance R s from the center of the body with a constant angular velocity 

O = (/j. /R s ^)'*‘ // ^, where / a is the universal gravitational constant multiplied by the body 
mass and where the magnitude of Rg is [R s (t) • Rg(t)] l/2 . Consider an inertial 
XYZ-axis system fixed in the center of the body such that, at the initial time t Q , it is 
alined with the rotating xyz-axis system. In this framework the target can be pictured 
as in figure A-2. 
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Figure A-2.- Target viewed in the inertial axis system. 


The angles t Q and 6 0 define the normal and line of nodes, respectively, of the 
target orbital plane relative to the inertial system. The x' and y’ axes therefore 
define the orbital plane of the target. If at t 0 the target is in the position 
(x',y') = (R s cos cp Q , R s sin <p 0 ) and moves toward the line of nodes, then 


R s [x’(t),y’(t),z’(t)] 



COS ^Pq - 0(t - to)] 

sin(^ 0 - f2(t - t 0 )] > 


(A2) 


and 


R S (X,Y,Z) = TjRsCx’.ySz’) 


(A3) 


where 


T l = 


cos 9 0 
sin e Q 
0 


-cos t 0 sin 6 0 
cos t Q cos $o 
sin L 0 


sin t 0 sin 6 0 
-sin to cos #o 
cos l q 


(A4) 


Since the xyz-axis system rotates about the Z-axis with a constant angular 
velocity co, 

R s (x,y,z) = R s (t) = T 2 (t)R s (X,Y,Z) (A5) 
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where 


Also 


T 2 (t) = 


cos w(t - to) 
-sin co(t - to) 


0 


sin o)(t - to) 0 

cos co(t — t Q ) 0 

0 1 


dR s (t) 

dt 


= T 2 (t)R s (X,Y,Z) + T 2 (t)R s (X,Y,Z) 


(A6) 


where 


and 


T 2 (t) = 


co 


-sin co(t - t Q ) 
-cos co(t - t Q ) 


0 


cos co(t - t Q ) 0 

-sin co(t - t Q ) 0 

0 0 


sin |^o " - to)] 

R S (X,Y,Z) = RsOTW-cos^o - n(t - t Q )]| 


Thus the position and rate of R s (t) can be obtained by specifying R s , l q , 6 0 , 
and cp Q at t 0 and by using equations (A5) and (A6). 

The thrust control vector T is related to the rotating axis system by 



cos 6q cos (pc' 

T = T| cos 6 C sin <p c 

sin 9 C 

as shown in figure A- 3. 


(A7) 


The vectors i, j , and k are unit vectors 
in the direction of the x, y, and z axes, respec- 
tively, and T is the magnitude of the thrust vector. 
Let 



~ u f 


cos 6 C 

cos <Pc 

u = 

u 2 


cos 6 C 

sin cp c 


_ u 3_ 


sin O c J 
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u 4 = T 


and 


r = 


r x 


R v x “ R s x 


RVy " R 

V v z 


s y 

S Z 


r = 


A x 

*y 





X 1 = r x 

(relative 

x distance) 

x 2 = r x 

(relative 

x velocity) 

x 3 = r y 

(relative 

y distance) 

li 

X 

(relative 

y velocity) 

x 5 = r z 

(relative 

z distance) 

x 6 = *z 

(relative 

z velocity) 

Xrj = m(t) 

(instantaneous vehicle 

x g = t 

(time) 



In this framework the dynamic equations can be written as (ref. 5) 

(x(t 0 ) = x 0 ; x(tf) = o) 


u.Mu _ _ 

X = ^ +Y ( X ’ X 8) 


x rj = - (x 7 (t 0 ) = m(t 0 ) = m 0 ) ) (A8) 


x 8 = 1 


(x 8 (t 0 ) = t 0 ) 


where 


with 


\ _ _ 

STRe 0 r 

s fl' 


Ax + R s || 2 ^- 


+ £2 2 MR s (x 8 ) + (n’ + 2wK + co 2 l)x (A9) 


x = COl^Xj, . . . x g ) 
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In order to compute the initial value of x(t 0 ), the initial value of R v (t 0 ) 
specified by 


can be 


Ry (to) — 


cos e v ° 

cos cp v ° 

cos e v ° 

sin cp v ° 

sin 6 > v ° 


where 0 V ° and cp y ° are as shown in figure A-4. Also 


z 



Figure A-4.- Initial orientation of vehicle with 
respect to the rotating axis system. 


Ry (to) 


Rv x (to) 

R-Vy (to) 

R Vz (to) 


and R s (to) and R s (to) are computed from 
equations (A2) to (A6). 


The act of rendezvous requires that the vehi- 
cle and target have the same position and velocity 
at tf ; hence, the condition x(tf) = 0. In addition, 
U4 = /3, the largest value obtainable for the thrust 
magnitude. 
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PROGRAM LISTING 

The program presented on the following pages is written in FORTRAN IV language 
for the Control Data series 6000 computer systems at the Langley Research Center. 
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PROGRAM E 1 257 ( I NPUT . OUTPUT « TAPES* I NPUT * TAPE6=OUTPUT ) 
C 

C TIME OPTIMAL RENDEZVOUS STUDY 

C 
C 

C DEFINITIONS 

C 

C NO 

C SOMEG 

C BETA 

C C 

C TF 

C PHIVO.THETVO 


CASE NO. 

ANG VEL OF CENTRAL BODY ABOUT ITS OWN AXIS 
MAX. THRUST 

EFFECTIVE EXHAUST SPEED 
FINAL TIME - GUESS 
INITIAL POSITION OF VEHICLE 


C RV 

C DRVXO * DRV YO « DRVZO 

C PH 1 0 • THET A 0 ♦ I 0 

C RS 


MAG. 

RVXO 

INIT 

MAG. 


OF RAD. VECTOR TO VEHICLE 
DOT. RVYO DOT. RVZO DOT 
AL POSITION OF ORBITING VEHICLE 
OF RAD. VECTOR TO ORBITING VEHICLE 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


VAR ( 1 ) 

VAR ( 8 ) 

MU 

VAR ( 9 ) —VAR ( 15) 
Cl 

SPEC 


I PR I NT 

IEROR 

IMAT 

LAMBDA 

CRIT 

MAXIT 

B 

INPUT AS FOLLOWS 
INPUT CARD NO. 


1 

2 

3 

4 

5 

6 

7 

8 
9 


INITIAL 

INITIAL 


I 1 Mfc 
MASS 


GRAVITATIONAL CONSTANT OF THE CENTRAL BODY 
I .C. ON PSll - PS 17 
COMPUTING INTERVAL 

SPEC = 1 • OE I 0 DO NOT PRINT AT SUB— I NTERV ALS 

SPEC .LT. 1.0E10 USE AS PRINT INTERVAL 

1 * PRINT PART I ALS - 0 * NO PRINT 

0*DO NOT PRINT TRUNCATION ERRORS - 1=PRINT 

1 * PRINT MATRICIES - 0 * NO PRINT 

PARAMETER USED IN ITERATION SCHEME 

STOPPING CRITERION 

MAX NO. OF ITERATIONS DESIRED 

DIAGONAL MATRIX OF WEIGHTING FACTORS 


QUANTITY FORMAT 

NO I 

SOMEG « BET A ♦ C . TF E 

PHI VO.THETVO.RV E 

DRVXO ♦ DR VYO . DRVZO E 

PH 10. THET AO ♦ IO.RS E 

VAR ( 1 ) . V AR < 8 ) * MU E 

VAR < 9 ) - VAR (12) E 

VAR (13) - VAR (15) E 

Cl. SPEC E 
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c 


10 


IPRINT. IEROR. IMAT 

I 


c 


1 1 


LAMBDA • CR I T * MAX I T 

E • I 


c 


12 


B ( 1 ) -B ( 4 ) 


E 


c 

c 


13 


B<S)-B<7) 


E 


c 

COMMON 

/SPACE/ 







1 

VAR 

•CUVAR 

• SAVE 

♦ C 

• MAX IT 

• I mat 


2 

E 

♦ PE 

•PEMAT 

•PEVEC 

•ERRVAL 

• SOMEG 


3 

DER 

• RS 

• BETA 

• DPI 

• TF 

• ELE1 


4 

RSV 

• FI 

• DP2 

• Cl 

• ELE2 

• DRSV 


5 

F2 

• El 

• I I 

♦ SOMEGS 

•COMEGS 

• IPRINT 


6 

F3 

• E2 

• N 

• COMEG 

• EN2 

•KOUNT 


7 

Q1 

• E3 

•tempci 

• SIO 

• PHIO 

• Q2 


a 

E4 

• SPEC 

• STO 

•LAMBDA 

•TEMPT 

• Q3 


9 

T2 

•TEMPSP 

• CIO 

• CRIT 

♦ EN 1 

♦ I KOUNT 


1 

RVXO 

• RVYO 

♦ RVZO 

• DRVXO 

♦DRVYO 

• DRVZO 


2 

T4 

• MU 

• CTO 

♦ 0 

♦ XT 1 



C 

REAL LAMBDA* MU 
REAL ! 0 
LOGICAL FIRST 
C 
C 


D I MENS I ON 

1 VAR ( 93 ) 

2 ELE1 (92 ) 

3 RSVC3) 

4 PPSI (6.6 ) 

5 CUPS 1(6*6) 

6 DRPS I (6*6) 

7 ERPS 1(6*6) 

8 PE (7*7) 

9 B ( 7 ) 

C 

C 


•CUVAR ( 93 ) 
*ELE2 (92 ) 

• DRSV (3 ) 

« PPS 17(6 ) 

♦ CUPS 17 ( 6 > 
♦DRPS 17(6) 
•ERPS 17(6) 

* PEMAT (7*7) 
• SAVE ( 93 ) 


* DER ( 93 ) 
•ERRVAL (92 ) 

* PX (6*6) 

* CUPX (6*6 ) 

* DRPX ( 6*6 ) 

♦ ERX (6*6) 

* E ( 7 ) 

♦PEVEC (7.1) 
*XT1 ( 3 * 3 ) 


EQUIVALENCE 

1 (VAR ( 16) *PX( 1*1)) 

( VAR ( 88 ) * PPS 17(1 ) ) 

( CUVAR ( 52 ) ♦ CUPS 1(1*1)) 

( DER (16) *DRPX ( 1*1)) 

(DER ( 88 ) *DRPS I 7 ( 1 ) ) 
(ERRVAL (51 ) *ERPS I < 1.1)) 


« ( VAR ( 52 ) «PPSI ( 1 « 1 ) > 

* ( CUVAR (16). CUPX (1*1)) 

• ( CUVAR ( 88 ) . CUPS 17(1)) 

. (DER (52 ) *DRPS I (1*1)) 

« ( ERRVAL (15). ERX (1.1)) 

. (ERRVAL (87), ERPSI7(1 )) 


C 


EXTERNAL DERSUB.CHSUB 
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EQU I VALENCE ( VAR ( 1 > • T 1 ) 

FORMAT STATEMENTS 

1 00 FORMAT ( I 20/4E20 .8/3E20 • 8/3E20 . 8/4E20 . 8/3E20 • 8/4E20 • 8/3E20 • 8 ) 

200 FORMAT < 2E20.8/3 1 20/2E20 .8.1 20/4E20 .8/3E20 .8 > 

500 FORMAT ( 30H1 T I ME OPTIMAL RENDEZVOUS STUDY. XOX 8HCASE NO. 13////// 
16H INPUT///5H BETA 2X E16.8.10X 2HRS 8X E16.8, 1 OX 7H0MEGA M 3X 

2E16.8. 1 OX 6HLAMBDA 2X E16.S/2H C 5X E16.8. 1 OX 2HRV 8X E16.8. 1 OX 

37H0MEGA T 3X E16.8. 1 OX 2HMU 6X E16.8//) 

50! FORMAT ( 1 9H WEIGHTING ELEMENTS/7E1 6.7/// ) 

502 FORMAT < I 9H INITIAL CONDI T I ONS//3H TO .4X . El 6 . 8 . 1 OX .6HPH I VO. AX. 

1 El 6.8. 1 OX . 8HRVX0 DOT , 2X . El 6 .8. 1 OX » 4HPH 10. 4X « E 1 6 . 8/5H MASS . 2X . El 6 • 8 
2.1 OX « 8HTHET A VO . 2X . El 6 . 8 1 1 OX » 8HRVYO DOT . 2X . El 6 .8 » 1 OX . 6HTHETA0 . 2X . 
3E16.8/5H MD0T2XE1 6.8, 46 X 8HRVZ0 D0T2X E 1 6 .8 . 1 0X2H I 0 6XE16.8///) 

503 FORMAT ( 1 9H ASSUMED COND I T I ONS//3H TF.4X.E16. 8* 10X . 4HPS I 1 .6X.E16.8. 
1 10X.4HPSI2.6X.E16.8/5H PS 1 7 . 2X . El 6 . 8 » 1 OX . 4HPS 1 3 . 6X . E 1 6 . 8 . 1 OX. 
24HPSI4.6X.E16.8/33X.4HPSI5.6X.E16.8. 10X.4HPS I6.6X.E16.8/// ) 

504 FORMAT (28H COMPUTED INITIAL COND I T I ONS//3H XR . 4X . El 6 .8 . 1 OX . 6HXR DO 
1T.4X.E16.8/3H YR . 4X .E 1 6 .8 . 1 OX . 6HYR DOT . 4X ,E 1 6 . 8/3H ZR.4X.E16.8. 
210X.6HZR DOT * 4X ♦ E 1 6.8//4H RSX . 3X.E1 6.8 « 1 OX .7HRSX DOT . 3X . El 6 .8 , 1 OX . 
33HRVX.7X.E16.8/4H RSY . 3X. El 6. 8 ♦ 1 OX ♦ 7HRSY DOT , 3X . El 6. 8 . 10X.3HRVY. 
47X.E16.8/4H RSZ . 3X . El 6 • 8 . 1 OX . 7HRSZ DOT . 3X . E 1 6 .8 ♦ 1 OX . 3HRVZ . 7X . E 1 6 . 8 
5/53H1 INTEGRATION ROUTINE USES FIXED COMPUTING INTERVAL = , E16.8//) 

800 FORMAT { 5H TIME 2X E16.8. 1 OX 4HRELD 5X El 6.8, 1 0X6HTHRUST 3X E16.8/ 
15H MASS 2X E16.8. 10X 4HRELV 5X E 1 6. 8 « 1 0X2HRV 7X E16.8/) 

804 FORMAT (5H PS 1 1 2X E16.8.10X 4HPSI2 5X E16.8.10X 4HPSI7 5X E16.8/ 

15H PS 1 3 2X El 6.8. 1 OX 4HPSI4 5X E16.8/5H PS 1 5 2X E16.8. 10X 4HPSI6 

25X El 6. 8/ ) 

805 FORMAT (9H PART I ALS//8H X/ALPHA/6 ( 6E 1 8 . 8/ ) // 

1 1 OH PS I /ALPHA/6 (6E1 8. 8/ ) //I 1 H PS 1 7/ALPHA/6E 1 8 . 8// ) 

806 FORMAT { 3H XR . 4X . E 1 6 .8 * 1 OX . 6HXR DOT ♦ 3X , E 1 6. 8 . 1 OX « 7HDXR DOT.2X. 
1E16.8. 10X.2HUX.8X.E16.8/3H YR.4X.E16.8. 10X.6HYR DOT . 3X . E 1 6 . 8. 1 OX . 
27HDYR D0T.2X.E16.8. 10X.2HUY.8X.E16.8/3H ZR . 4X . El 6 .8 , 1 OX « 6HZR DOT. 
33X . E 1 6 . 8 . 1 OX . 7HDZR DOT . 2X . E 1 6 . 8 . 1 OX . 2HUZ . 8X . E 1 6.8/ ) 

808 FORMAT <4H RSX . 3X . E 1 6 . 8 , 1 OX . 7HRSX DOT . 2X . E 1 6 . 8 . 1 OX . 3HRVX . 6X . E 1 6 . 8 . 

1 1 OX . 7HRVX DOT » 3X . E 1 6. 8/4 H RSY . 3X. E 1 6 . 8 . 1 OX ♦ 7HRSY DOT , 2X . E 1 6 .8 , 1 OX . 
23HRVY.6X.E16.8. 10X.7HRVY DOT , 3X . E 1 6 . 8/4H RSZ . 3X . El 6 . 8 . 1 OX ♦ 

37HRSZ D0T.2X.E16.8, 1 OX » 3HRVZ ♦ 6X . E 1 6 . 8 , 10X.7HRVZ DOT . 3X . E 1 6 . 8/ ) 

810 FORMAT { 5H TIME 2X E16.8/1 

811 FORMAT ( 5H PS 1 1 2X E16.8.10X 4HPSI2 5X E16.8.10X 4HPSI7 5X E16.8. 

UOX 2HUX 8X E16.8/5H PSI3 2X E16.8. 1 OX 4HPSI4 5X E16.8. 45X 

2 2 HUY 8X E16.8/5H PSI5 2X E16.8. 1 OX 4HPSI6 5X E16.8, 45X 2HUZ 8X 

3E16.8/) 
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900 FORMAT (24 H LOCAL TRUNCATION ERRORS//2H X/6E18.8//3H X7/E18.8// 

1 AH PSI/6E18.8//5H PS I 7/El 8. 8// II H P (X/ ALPHA )/6 (6E1 8. 8/ )// 

213H P CPS I /ALPHA )/6 < 6E 1 8 •8/ ) // 1 4H P C PS 1 7/ALPHA ) / 

36E18.8//) 

4900 FORMAT (//) 

901 FORMAT < 40H MAX NO* OF ITERATIONS HAS BEEN REACHED) 

903 FORMAT { 9H ECALPHA) 2X E16.8/) 

START 

READ INPUT 

1 READ(5« lOO ) NO » SOMEG. BETA »C « TF • PH I VO « THETV0« R V. DRVXO • DRVYO .DRVZO « 
1 PH 10* THET AO « IO.RS.VARC 1 ) .VAR(8 ) .MU. < VAR ( I ) . 1=9. 15) 

2 READ (5« 200 )C I « SPEC • I PR t NT . I EROR . I MAT « LAMBDA « CR I T .MAX I T « 

1 (BCI )« 1=1 *7) 

COMPUTE CONSTANT TERMS AND INITIAL CONDITIONS 

COMEGS*MU/RS**3 
COMEG* SORT ( COMEGS ) 

SOMEGS*SOMEG**2 
SIO*SINC 10 ) 

STO*S I N (THET AO ) 

CIO*COS< 10 ) 

CTO*COS CTHETAO ) 

T 1 MATRIX 

XT 1 (1.1 ) =CTO 
XT 1 Cl «2>=-Cl0*ST0 
XT 1 C 1 . 3 ) =S I 0*ST0 
XT! (2. 1 ) =STO 
XT l (2.2)=CI0*CT0 
XT 1 (2.3)=-SI0*CT0 
XT 1 (3*1 )*0*0 
XT1 (3.2 >=SI0 
XT 1 (3.3 >=CIO 

INITIALIZATION 

RVX*RV*COS ( THET VO ) *COS ( PHI VO ) 

RVY*RV*COS (THET VO ) *S I N (PHI VO ) 

RVZ*RV*S I N ( THET VO ) 
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RVXO*RVX 
RVYO*RVY 
RVZO*RVZ 
CALL COWP(O.O) 

VAR (2 ) *RVX-RSV ( 1 ) 

VAR (3 )*DRVXO— DRSV ( 1 ) 

VAR (4 )*RVY-RSV ( 2 ) 

VAR ( 5 )*DRVYO-DRSV ( 2 ) 

VAR ( 6 ) sRVZ-RSV ( 3 ) 

VAR ( 7 ) -ORVZO-DRSV < 3 ) 

DER <8 )— BETA/C 

WR I TE ( 6 » 50 0 ) NO ♦ BETA * RS • SOMEG * LAMBDA , C » R V , COMEG * MU 
WRITE (6.501 ) <B(I).I*1*7) 

WRITE(6.502 ) VAR ( 1 ) .PHI V0.DRVX0.PHI0.VAR (8 ) . THETVO . DRVYO.THETAO, 

1 DER ( 8 ) • DRVZO .10 

WRITE (6. 503 ) TF ♦ V AR ( 9 ) « VAR (10). VAR ( 1 5 ) « VAR ( 1 1 ) « VAR (12). VAR (13), 

1 VAR ( 14 ) 

WRITE (6, 504) (VAR ( I ). 1*2.7) ,RSV( 1 ) ,DRSV( 1 ) .RVX.RSV(2 ) .DRSV(2 ) ,RVY, 
2RSV ( 3 ) . DRSV ( 3 ) , RVZ , C I 
C 

C INITIALIZATION FOR INTEGRATION ROUTINE 

C 

C 

DO 17 1*1.15 
17 SAVE ( I )*VAR( I ) 

DO 14 1*16.93 
14 VAR ( I )*SAVE( I ) 

C 

TEMPC I =C I 
TEMPSP*SP£C 
I KOUNT*0 
1010 KOUNT=0 
11*0 

FIRST* .TRUE. 

CI*TEMPCI 
SPEC*TEMPSP 

10 CALL INT2C 1 1 .N.NT « C I .SPEC.C I MAX, I ERR • VAR , CUVAR, DER.ELE1 , ELE2.ELT. 

1 ERR VAL « DERSUB , CHSUB « I TEXT ) 

RETURN FROM INTEGRATION ROUTINE 

11 IF ((ABS(Tl-TF) .LE. 1.0E-06) .OR. FIRST .OR. 

1 ( (SPEC .LT. 1.0E10) .AND. (SPEC .NE. 0.0))) GO TO 16 

GO TO 20 
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WRITE OUTPUT 

16 RELD=SQRT < VAR ( 2 )**2+VAR (4 )**2+VAR (6 )**2 > 

RELV*SQRT ( VAR ( 3 ) **2+VAR ( 5 ) **2+ VAR 1 7 ) **2 ) 

RVX= VAR ( 2 ) +RSV ( 1 > 

RVY= VAR ( 4 ) +RSV ( 2 ) 

RVZ* VAR C 6 ) +RSV ( 3 ) 

RVMAG=SORT<RVX**2+RVY**2+RVZ**2 ) 

0RVX*VAR(3 l+ORSVU ) 

DRVY = VAR ( 5 ) +DRS V ( 2 ) 

DRVZ*VAR < 7 l+DRSV ( 3 ) 

U4»BETA 
UX*VAR( 101/EN1 
UY=VAR( 12)/EN1 
UZ * VAR (14) /EN 1 
IF (FIRST) GO TO 60 

WRITE (6*800 ) VAR ( 1 ) .RELD.U4 * VAR (8 ) «RELV .RVMAG 

WR I TE ( 6 « 804 ) VAR ( 9 ) « VAR ( 1 0 ) ♦ VAR (15) ♦ VAR ( 11 ) «VAR(12)« 

1VAR(13) « VAR (14) 

WRITE(6«806) VAR ( 2 ) « VAR (3 ) . DER ( 3 ) ,UX« VAR ( 4 ) . VAR (5 > «DER (5 ) « UY . 
1 VAR { 6 ) ♦ VAR ( 7 ) « DER ( 7 ) « UZ 

WRITE(6«80S ) RSV ( 1 ) .DRSV( 1 ) «RVX «DRVX ,RSV ( 2 ) * DRSV (2 ) ,RVY«DRVY« 
1 RSV ( 3 ) « ORSV ( 3 ) .RVZ.DRVZ 
GO TO 61 

60 WRITE (6. 810) VAR ( 1 ) 

WR ITE ( 6 »8l 1 ) VAR (9 )♦ VAR ( 10 ) «VAR< 15) «UX« VAR( 1 1 ) »VAR( 12 ) ,UY* 

1 VAR (13). VAR ( 1 4 ) *UZ 

61 IF (IPRINT .NE. 1 ) GO TO SO 

WRITE (6. 80S) ( (PX (I.J).J=1.6)*1=1.6)« 

1 ( (PPSI ( I . J) . J*1 *6 ) . 1=1 .6) • ( PPS 17(1 ) « 1=1 «6 ) 

SO I F ( I EROR ,EQ. 1 )GO TO 9 
GO TO 49 

9 WRITE (6.900 ) (ERRVAL ( I ) . 1=1 . 14 ) . ( (ERX( I . J ) . J= 1 .6 ) . I = 1 .6 ) . 

1 ( (ERPSI ( I*J).J=1*6)«I=1 .6). (ERPS 17 ( I ) ♦ 1=1 .6 ) 

49 WRITE (6*4900 ) 

20 IF (FIRST) FIRST* .FALSE. 

IF (ABS(Tl-TF) .LE. 1.0E-06) GO TO 13 
IF ((Tl+CI) .LE. TF) GO TO 10 
C I *TF— T 1 
11=0 

SPEC=0*0 
GO TO 10 
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13 IKOUNT* IKOUNT+1 

I F ( I KOUNT .GT* MAXI T ) GO TO IS 

COMPUTE (E.8*E>/2 


C 

c 

c 


EDP*0.0 
DO 19 1*2*7 

19 EDP*£DP+B( 1-1 >*VAR( I)**2 

EOP* • 5* ( EDP+B < 7 ) * VAR ( 1 5 > **2 > 
WRITE <6 *903 > EDP 
IF (EDP .LE. CR1T ) GO TO 1 
CALL ITERAT 
DO 18 1 * 1.93 
18 VAR ( I ) *SAVE ( I ) 

GO TO 1010 
15 WRITE( 6.901 ) 

GO TO 1 
ENC 

SUBROUTINE DERSUB 
COMMON /SPACE/ 


I 

VAR 

♦CUVAR 

• SAVE 

* c 

♦ MAX IT 

, I MAT 

2 

E 

♦ PE 

•PEMAT 

•PEVEC 

♦ERRVAL 

♦ SOMEG 

3 

DER 

*RS 

♦ BETA 

♦ DPI 

♦ TF 

, ELE1 

4 

RSV 

• Fi 

♦ DP2 

• Cl 

♦ ELE2 

♦ DRSV 

5 

F2 

• El 

• 1 I 

•SOMEGS 

•COMEGS 

• IPRINT 

6 

F3 

• E2 

♦ N 

♦COMEG 

• EN2 

♦KOUNT 

7 

Q1 

• E3 

• TEMPCI 

♦ SIO 

• PH 10 

« Q2 

8 

E4 

♦ SPEC 

• STO 

♦LAMBDA 

•TEMPT 

.03 

9 

T2 

♦TEMPSP 

♦ CIO 

• CRIT 

• EN1 

, IKOUNT 

1 

RVXO 

♦ RVYO 

♦ RVZO 

• DR VXO 

•DRVYO 

, DRVZO 

2 

T 4 

• MU 

♦ CTO 

♦ B 

♦ XT 1 



REAL LAMBDA, MU 


D I MENS 1 ON 


1 

VAR (93) 

• CUVAR ( 93 ) 

, DER ( 93 ) 

2 

ELE1 (92) 

•ELE2 (92 ) 

.ERRVAL (92) 

3 

RSV ( 3 ) 

• DRSV <3 ) 

«PX( 6.6) 

4 

PPSI (6 ♦ 6 ) 

♦ PPS 17(6) 

. CUPX (6.6) 

5 

CUPSI ( 6 ♦ 6 ) 

♦ CUPS 1 7 C 6 > 

♦ DRPX (6.6) 

6 

DRPSI (6*6) 

♦DRPSI7C6) 

» ERX (6,6) 

7 

EPPS I (6 ♦6> 

♦ERPS 17(6) 

«E(7) 

8 

PE ( 7 « 7 ) 

♦ PEMAT ( 7^ 7 ) 

♦PEVEC (7. 1 ) 

9 

B ( 7 ) 

•SAVE (93) 

»XT1 (3.3) 
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EQUIVALENCE 

1 (VAR ( 16 > »PX( 1*1)) 

2 (VAR (88) ,PPSI7(1 >) 

3 ( CUVAR ( 52 ) * CUPS 1(1*1)) 

4 (DER(16) «DRPX( 1 « 1 ) ) 

5 ( DER ( 88 ) « DRPS I T ( 1 ) ) 

6 (ERRVAL (51 > *ERPSI ( 1 • 1 > ) 

C 

EQUIVALENCE (CUVAR ( 1 ) «TJ> 

c 

TCOMP“T J— SAVE ( 1 ) 

CALL COMP (TCOMP ) 

Q1=CUVAR(2)+RSV( 1 ) 

Q2“CUVAR ( 4 )+RSV ( 2 ) 

Q3=CUVAR (6 )+RSV (3 ) 

ENI =SQRT (CUVAR ( 10 )**2+ CUVAR ( 1 2 >**2+ CUVAR ( 14 ) **2 > 
EN2=SQRT(Q1**2+Q2**2+Q3**2 ) 

OP 1 “CUVAR ( 2 ) *CUVAR (10) +CUVAR ( 4 ) *CUVAR (12) +CUVAR ( 6 ) #CUVAR (14) 
DP2=RSV( 1 )*CUVAR ( 10)+RSV(2 )*CUVAR( I 2 )+RSV ( 3 ) *CUVAR ( 14 ) 

F2“ (COMEGS*RS**3)/EN2**3 
F3= (3«0*F2/EN2**2 )* (DPI +0P2 ) 

Et*l •0/’ ( CUVAR ( 8 ) *EN i ) 

E2“E1 /CUVAR ( 8 ) 

E3=E1/EN1**2 
E4“EN1 /CUVAR (8)**2 
T2* ( 3 • 0*F2 ) /EN2**2 
T4*5 • 0*F3/EN2**2 
C 

C STATE VARIABLES 

C 

C 

DER (9 )=F2*CUVAR ( 1 O ) — F3*Q 1 -SOMEGS*CUV AR (10) 

DER (10) “—CUVAR ( 9 > +2 • 0*SOMEG*CUVAR (12) 

DER (11 )“F2*CUVAR ( 1 2 ) -F3*Q2-SOMEGS*CUVAR ( 12 ) 

DER ( 1 2)“-2«0*SOMEG*CUVAR( 10 ) -CUVAR (11) 

DER (13 )“F2*CUVAR ( 14)-F3*Q3 
DER (14) --CUVAR (13) 

F 1 “BETA/ ( CUVAR ( 8 ) *EN 1 ) 

51 DER ( 2 > “CUVAR ( 3 ) 

DER ( 3 ) “F 1 *CUVAR (10) -F2*Q 1 +C OMEGS*RSV ( 1 ) +SOMEGS*CUV AR { 2 ) +2 . 0*SOMEG* 
1 CUVAR ( 5 ) 

DER ( 4 ) “CUVAR ( 5 ) 

DER (5 )“F1*CUVAR ( 12)-F2*Q2+COMEGS*RSV (2 ) -2« 0*SOMEG*CUVAR (31+SOMEGS* 


« (VAR (52 ) «PPSI ( 1 ♦ 1 ) ) 

« ( CUVAR ( 16) .CUPXl 1 t l ) ) 

» ( CUVAR ( 88 ) « CUPS 17(1 ) ) 

« (DER (52 ) » DRPS I ( 1 « 1 > ) 

* ( ERRVAL ( 1 5 ) « ERX ( 1 « 1 > ) 

« ( ERRVAL ( 87 ) tERPS 17(1 )) 
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II 


1 CUVAR (4 ) 

DER(6)»CUVAR(7) 

OER ( 7 ) »F1 *CUVAR (14) -F2*Q3+COMEGS*RSV ( 3 ) 
DER ( 8 ) *— BE T A /C 

DER(15)= (BETA*ENl )/CUVAR(8 ) **2 


PARTIALS 

DO 60 1=1 «6 
DRPX (1*1 ) =CUPX ( 2 « I ) 

C 

DRPX <2 • I )«BETA*( ( El -E3*CUVAR ( 10 >**2 >*CUPSI (2*1 )-E3*CUVAR ( 1 0 )* 

1 CUVAR ( 1 2 ) *CUPS 1 (4.1 ) -E3+CUVAR ( 1 0 >*CUVAR ( 14)*CUPSI (6.1))+ (T2*Ql**2- 
2F2+S0MEGS ) *CUPX ( 1 .1 ) +T2*Q1 *Q2*CUPX ( 3 . I >+2.0*SOMEG*CUPX (4. I )+T2*01* 
3Q3*CUPX(5. I ) 

C 

DRPX (3.1 ) *CUPX ( 4 « I) 

c 

DRPX (4. I )=BETA* (-E3+CUVAR ( 1 0 ) *CUVAR ( 1 2 ) *CUPS 1 (2.1 )+(El-E3* 

1 CUVAR (12) **2 )*CUPSI (4.1 )-E3*CUVAR( 12)*CUVAR ( 14 )*CUPSI (6.1) )+T2*Ql* 
2Q2*CUPX(1 . I )-2.0*S0MEG*CUPX(2. I )+ ( T2*Q2**2— F2+S0MEGS ) *CUPX (3 . I )+ 
3T2*Q2*Q3*CUPX (5.1) 

C 

DRPX (5. I )=CUPX(6. 1 ) 

C 

DRPX (6.1 )=BETA* ( -E3*CUVAR (10) *CUVAR (14) *CUPS I ( 2 « I ) -E3* 

1 CUVAR (12 )*CUVAR ( 14 )*CUPSI (4.1 )+ (El-E3*CUVAR ( 1 4 >**2 )*CUPS I (6 . I ) )+ 
2T2*Q1*Q3*CUPX( 1 . 1 >+T2*Q2*Q3*CUPX ( 3 . I )+ (T2*Q3**2-F2 )*CUPX (5.1) 

C 

DRPSI (1 « I )=(F2-T2*Ql**2-SOMEGS)*CUPSI (2. I >-T2*Ql*Q2*CUPSI (4.1 )-T2* 
1Q3*Q1*CUPSI (6. 1 ) + ( —2. 0*T2*CUVAR (10 ) *Q 1 -F3+T4*Q1**2 ) *CUPX ( 1 . I )+(-T2 
2* ( CUVAR (10 )*Q2+CUVAR ( 1 2 )*Q1 )+T4*Ql*Q2 )*CUPX(3 . I )+ (-T2* (CUVAR ( 1 0 ) * 
3Q3+CUVAR (14) *Q 1 )+T4*Q 1 *Q3 )*CUPX ( 5 . I ) 

C 

DRPSI (2.1 )=— CUPS! (1.1 )+2.0*SOMEG*CUPSI (4.1) 

C 

DRPSI (3.1 )=-T2*Qi*Q2*CUPSI (2. I )+ (F2-T2*Q2**2-SOMEGS )*CUPS1 (4. I )- 
1T2*Q2*Q3*CUPSI (6.1 ) + (-T2* (CUVAR ( 12 )*Q 1 +CUVAR ( 1 0 >*Q2 > + T4*Q1 *Q2 ) * 
2CUPX ( 1 « I )+(-2.0*T2*CUVAR( 12 )*02-F3+T4*Q2**2 )*CUPX(3. I >+ (-T2* ( 

3 CUVAR (12) 4Q3+CUVAR (14) *Q2 ) +T4*Q2*Q3 ) *CUPX (5.1) 

C 

ORPSI (4.1 >=-2.0*SOMEG*CUPSI (2.1 ) -CUPS I (3.1) 

C 

DRPSI (5.1 )=-T2*Ql*Q3*CUPSI (2,1 )-T2*Q2*Q3*CUPSl (4.1 )+(F2-T2*Q3**2 )* 
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1 CUPS I (6*1 >+ ( — T2* ( CUVAR < 1 4 ) *Q1 +CUVAR (10 >*Q3 )+T4*Ql*Q3 )*CUPX ( 1 * I ) + 
2 ( -T2* ( CUVAR (14) *Q2+CUVAR (12) *Q3 ) +T4*Q2*Q3 ) *CUPX ( 3 ♦ I ) + ( -2 * 0*T2* 
3CUVAR (14) *Q3-F3+T4*Q3**2 ) *CUPX (5*1) 

C 

DRPSI (6* ! )=— CUPS1 (5* I ) 

C 

DRPS I 7 ( I ) =BETA*E2* ( CUVAR ( 10) *CUPSI (2*1 >+CUVAR( 12 )*CUPSI (4*1) 

1 +CUVAR (14) *CUPS I (6*1)) 

C 

60 CONTINUE 
RETURN 
END 

SUBROUTINE CHSUB 
COMMON /SPACE/ 


1 

VAR 

# CUVAR 

• SAVE 

♦ C 

♦MAX IT 

♦ I MAT 

2 

E 

• PE 

•PEMAT 

•PEVEC 

tERRVAL 

« SOMEG 

3 

OER 

*RS 

• BETA 

♦ DPI 

*TF 

• ELE1 

4 

RSV 

• FI 

♦ DP2 

*CI 

* ELE2 

• DRSV 

5 

F2 

♦ El 

♦ I I 

•SOMEGS 

* COMEGS 

« I PR I NT 

6 

F3 

♦ E2 

♦ N 

•COMEG 

♦ EN2 

.(COUNT 

7 

01 

♦ E3 

♦ TEMPO I 

♦ SIO 

• PHI 0 

• 02 

8 

E4 

♦ SPEC 

♦ STO 

•LAMBDA 

•TEMPT 

«Q3 

9 

T2 

•TEMPSP 

• Cl 0 

• CR1T 

• EN1 

• IKOUNT 

1 

RVXO 

• RVYO 

♦ RVZO 

♦DRVXO 

•DRVYO 

* DRVZO 

2 

T 4 

♦ MU 

• CTO 

♦ B 

• XT 1 


REAL 

LAMBDA# MU 







C 

C 


DIMENSION 



1 

VAR < 93 ) 

♦ CUVAR ( 93 ) 

• DER ( 93 ) 

2 

ELE1 (92 ) 

• ELE2 (92 > 

•ERRVAL (92 ) 

3 

RSV ( 3 ) 

• DRSV ( 3 > 

♦PX (6*6) 

4 

PPSI (6 • 6 ) 

•PPS 17(6) 

* CUPX (6*6) 

5 

CUPSI < 6 «6 ) 

* CUPS 17(6) 

• DRPX (6*6 > 

6 

DRPSI ( 6 ♦S ) 

•DRPS I 7 ( 6 ) 

♦ERX(6*6 ) 

7 

ERPSI (6*6) 

. ERPS 17(6) 

• E(7 ) 

8 

PE ( 7 • 7 ) 

. PEMAT (7.7) 

♦PEVEC (7 • 1 ) 

9 

B ( 7 ) 

•SAVE (93) 

♦ XT 1 ( 3 • 3 ) 


c 

c 

equivalence 

1 (VAR (16) «PX (1*1)) * (VAR (52 ) .PPSI ( 1 ♦ 1 > ) * 

2 (VAR (88) »PPSI7( 1 ) ) ♦ ( CUVAR ( 1 6 ) .CUPX ( 1 « 1 ) ) ♦ 

3 ( CUVAR ( 52 ) * CUPS 1(1*1)) .( CUVAR ( 88 ) .CUPSI7 ( 1 ) ) * 
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4 (DER(16> *DRPX( 1 * 1) ) 

5 (DER(88) *DRPSI7( 1 > > 

6 (ERRVAL (51 ) ♦ ERPS 1(1*1)) 

EQUIVALENCE (CUVAR ( 1 ) *T I ) 


* ( DER ( 52 ) «DRPS 1(1*1 ) ) « 
♦ (ERRVAL (15) *ERX(1 * 1 ) ) « 
♦(ERRVAL(87)*ERPSI7(1 ) ) 


IF ((TI+Ct) *LE. TF) GO TO 78 
11*3 

78 RETURN 
END 

SUBROUTINE COMP(DT) 

COMMON /SPACE/ 


1 

VAR 

•CUVAR 

♦ SAVE 

♦ C 

♦MAXIT 

♦ I MAT 

2 

E 

• PE 

♦PEMAT 

♦PEVEC 

♦ERRVAL 

* SOMEG 

3 

DER 

• RS 

♦ BETA 

♦ DPI 

♦ TF 

♦ ELE1 

4 

RSV 

• FI 

♦ DP2 

♦ Cl 

♦ ELE2 

♦ DRSV 

5 

F2 

♦ El 

♦ I I 

♦SOMEGS 

♦COMEGS 

♦ IPRINT 

6 

F3 

• E2 

♦ N 

♦COMEG 

• EN2 

♦KOUNT 

7 

QX 

• E3 

♦ TEMPC I 

• SIO 

♦ PH 10 

♦ Q2 

8 

E4 

• SPEC 

♦ STO 

♦LAMBDA 

•TEMPT 

♦ Q3 

9 

T2 

♦TEMPSP 

♦ CIO 

• CRIT 

• EN1 

« KOUNT 

1 

RVXO 

• RVYO 

♦ RVZO 

•DRVXO 

•DRVYO 

♦ DRVZO 

2 

T 4 

• MU 

♦ CTO 

♦ B 

♦ XT 1 



REAL LAMBDA* MU 


1 

D I MENS I ON 

VAR (93) 

♦ CUVAR ( 93 ) 

♦ DER ( 93 ) 

2 

ELE1 (92) 

♦ELE2<92) 

♦ERRVAL (92) 

3 

RSV (3) 

• DRSV ( 3 ) 

• PX ( 6 «6 ) 

4 

PPSI (6»6) 

♦PPS 17(6) 

♦ CUPX ( 6 «6 ) 

5 

CUPS I (6*6) 

♦CUPSI7(6> 

♦DRPX ( 6 ♦ 6 ) 

6 

DRPS 1 ( 6 ^6 > 

♦ DRPS 17(6) 

♦ ERX (6^6) 

7 

ERPSI (6 *6) 

♦ERPS 17(6) 

♦ E ( 7 ) 

8 

PE < 7 ♦ 7 ) 

♦ PEMAT ( 7« 7 ) 

•PEVEC (7^1 ) 

9 

B ( 7 ) 

•SAVE (93) 

* XT 1 < 3 ♦ 3 ) 


EQUIVALENCE 

1 (VAR (16) *PX< 1 « 1 > ) 

2 (VAR ( 88 ) «PPS 17(1 ) ) 

3 ( CUVAR ( 52 ) ♦ CUPS 1(1*1 ) > 

4 ( DER (16)* DRPX ( 1 * 1 ) > 

5 ( DER ( 88 ) » DRPS 17(1)) 


♦ (VAR (52 ) *PPSI (1*1)) 

« ( CUVAR (16) *CUPX ( 1 * 1 ) ) 
« ( CUVAR ( 88 ) ♦ CUPS 17(1 ) ) 
♦ ( DER ( 52 ) « DRPS 1 (1*1)1 
, (ERRVAL (15 ).ERX(1*1 ) ) 
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6 < ERRVAL. ( 5 1 ) « ERPS I < 1 « 1 ) > , < ERRVAL. ( 87 ) .ERPSI 7 < 1 > > 

DIMENSION XTO ( 3 ) « XT2 ( 3»3)«XT3(3) *XTD2<3»3 > .RSXYZ (3 ) »TDRS<3 ) * 
1 RSXYZD C 3 ) « TRSDOT < 3 > 

PH I OMT=PH I 0-COMEG*DT 

TO MATRIX 

XTO(l )=RS*COS (PHIOMT ) 

XTO I 2 ) »RS*S INI PHI OMT ) 

XT0(3>*0.0 

T2 MATRIX 

SOT*S I N ( SOMEG*DT ) 

COT»COS <S0MEG*DT ) 

XT2 ( 1 « 1 >®COT 
XT2 ( 1 »2>*S0T 
XT2 (1.3) *0.0 
XT2 (2. 1 ) =— SOT 
XT2 (2*2 ) *C0T 
XT2(2.3 >*0.0 
XT2I3. 1 )*0»0 
XT2 (3*2) =0 »0 
XT2(3.3)=1 .0 

T2 DOT MATRIX 

XTD2 (1.1) »-SOMEG*SOT 
XTD2 ( 1 »2)-S0MEG*C0T 
XTD2 ( 1 .3 >*0.0 
XTD2 (2.1) *-SOMEG*COT 
XTD2 (2.2) *— SOMEG*SOT 
XT D2 (2«3)=0»0 
XTD2 ( 3 « 1 )*0.0 
XTD2<3.2>*0.0 
XTC2<3.3>*0.0 

T3 MATRIX 

XT3 ( 1 >*C0MEG*XT0<2> 

XT3 (2 )«— COMEG* XTO ( 1 ) 

XT3 (3 ) *0.0 

TITO MATRIX 
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DO 100 1*1*3 
RSXYZ( 11*0*0 
DO 100 J*1 *3 

100 RSXYZt 1 ) *RSXYZ ( I )+XTl ( I »J)*XTO(J) 
RSV MATRIX 


DO 101 1*1*3 

RSV ( I 1*0*0 
DO 101 J*1 « 3 

1 01 RSV ( 1 ) *RSV ( I 1+XT2 ( I * J ) *RSXYZ ( J 1 
DO 102 1*1*3 

TORS ( I 1=0*0 
DO 102 U= 1*3 

102 TDRSt I )*TDRS (I >+XTD2( I . J)*RSXYZ< J) 

DO 103 1*1*3 

RSXYZDl 1 1=0*0 
DO 103 J=1 *3 

103 RSXYZDt I )=RSXYZD( I l+XTl ( I • J)*XT3(J) 

DO 104 1*1.3 

TRSDOT ( 1 1=0*0 
DO 104 J*1 *3 

1 04 TRSDOT ( 1 1 *TRSDOT ( I 1+XT2 < I « J 1 *RSXYZO ( J 1 


RSV DOT MATRIX 


DO 105 1=1*3 

105 DRSV (II =TDRS < I ) +TRSDOT (II 
RETURN 
END 

SUBROUTINE I TERAT 
COMMON /SPACE/ 


1 

VAR 

♦CUVAR 

• SAVE 

• C 

♦MAXIT 

• I MAT 

2 

E 

• PE 

•PEMAT 

•PEVEC 

«£RRVAL 

• someg 

3 

OER 

• RS 

• BETA 

• DPI 

» TF 

• ELE 1 

4 

RSV 

• FI 

♦ DP2 

• Cl 

♦ ELE2 

• DRSV 

5 

F2 

• El 

• 1 I 

•SOMEGS 

•COMEGS 

• IPRINT 

6 

F3 

♦ E2 

♦ N 

•COMEG 

* EN2 

•<OUNT 

7 

01 

♦ E3 

•TEMPCI 

• SIO 

*PH I 0 

• 02 

8 

E4 

• SPEC 

• STO 

•lamboa 

.TEMPT 

• Q3 

9 

T2 

•TEMPSP 

• CIO 

• CRIT 

«ENl 

♦ IKOUNT 

1 

RVXO 

• RVYO 

• RVZO 

•DRVXO 

*DRVY0 

• DRVZO 

2 

T 4 

• MU 

• CTO 

• B 

♦ XT 1 



C 

real lambda* mu 
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C 

c 


c 

c 


c 

c 

c 

c 

c 

c 


c 

c 


c 


DIMENSION 

1 VAR { 93 ) 

2 ELEl (92) 

3 RSV t 3 ) 

4 PPS I (6*6) 

5 CUPS I (6*6) 

6 DRPS I (6*6) 

7 ERPS 1(6*6) 

8 PE C7»7 ) 

9 B (7 ) 


t CUVAR ( 93 ) 

♦ ELE2 (92 > 
♦DRSV ( 3 ) 
♦PPSI7 (6 ) 

» CUPS 17(6) 
•DRPS 17(6 > 
*ERPSI7 (6 > 
♦PEMAT (7*7 ) 
* SAVE (93) 


EQUIVALENCE 

1 (VAR (16) *PX( 1 *1 > > • 

2 (VAR ( 88 ) «PPS 17(1 ) ) . 

3 ( CUVAR ( 52 ) * CUPS 1(1*1 > ) . 

4 (DER( 16 > »DRPX< 1 ♦ I > > * 

5 ( DER ( 88 > • DRPS I 7 ( 1 > ) ♦ 

6 (ERRVAL (51 ) *ERPS 1(1*1 1 > * 

DIMENSION I P I VOT ( 7 ) ♦ I NDEX (7*2) 


DIMENSION S AVMAT (7*7)* UNIT (7*7) 

FORM MATRIX OF PART I ALS 

DO 25 1*1 *6 
DO 25 J*1 »6 
PEC I*J)*PX(I .J) 

25 CONTINUE 


J *1 

DO 26 1*1 *5*2 
J=J+2 

PE ( 1 *7>*VAR( J) 

PE ( 1+1 «7)*DER( J ) 

26 CONTINUE 

DO 27 1*1 *6 

27 PE(7. 1 ) *PPS I 7 ( I > 

PE (7*7) *OER (15) 


« DER ( 93 ) 
♦ERRVAL (92) 
♦PX ( 6«6 ) 

♦ CUPX ( 6 ♦ 6 ) 
♦DRPX (6^6 ) 

♦ ERXC 6^6 ) 

♦ E <7 ) 

♦PEVEC(74 l > 

♦ XT1 ( 3#3 ) 


C VAR ( 52 )4PPS I ( 1 ♦! > ) 

( CUVAR ( 16) .CUPX( 1 ♦ 1 ) ) 

( CUVAR ( 88 ) ♦ CUPS 17(1 ) ) 
(DER (52) .DRPS I (1*1 )> 
(ERRVAL ( 15 ) tERX ( 1 * 1 ) ) 
(ERRVAL (87) *ERPS17(1 ) ) 



you 
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FORM E VECTOR 

DO 38 1 = 1.6 
38 E( I >=VAR( 1+1 ) 

E(7)*VAR(l5) 

C T 

C FORM PE *B*PE+LAM8DA*I MATRIX 

C 

C 

DO 40 1=1.7 
DO 40 J=1 .7 
PEMAT ( I . J)=0.0 
DO 40 K*1 .7 

PEMAT ( I « J >=PEMAT< I * U)+PE(K. I )*B (K >*PE(K. J ) 

40 CONTINUE 

DO 41 1*1 .7 

41 PEMAT ( I ♦ I >=PEMAT ( I . I l+LAMBDA 

C 

C T T 

C SOLVE THE MATRIX EQ (PE *B*PE+LAMBDA* I )DELTA=PE *B*E 

C 

DO 50 1*1.7 
PEVEC ( I . 1 >=0.0 
DO 50 J=1 .7 

50 PEVEC ( I . 1 )=PEVEC< I . 1 )-PE(U« I )*B(J)*E( J) 

IF ( 1 MAT .NE. 1) GO TO 13 
DO 14 1=1.7 
DO 14 J=1 .7 

14 SAVMAT ( I « J ) = PEMAT ( I « J } 

I F ( I MAT .EQ. 1) WR I TE (6.11 ) ( (PEMAT ( I . J > . J=1 .7 > , 1=1 .7 > 

11 FORMAT ( 6H PEMAT/7 ( 7E 16. 7/ )// > 

13 CALL MATINV(PEMAT, 7. PEVEC. 1 .DETERM. IPIVOT. INDEX.7. I SCALE > 
C 

IF ( I MAT .EQ. 1) WR I TE (6.12) ( ( PEMAT (I.J).J=1«7).I=1«7) 

12 FORMAT ( 8H INVERSE/7 (7E1 6.7/ )// ) 

IF ( I MAT .NE. 1> GO TO 15 

DO 16 1=1.7 
00 16 J=1 *7 
UNIT (I . J ) =0.0 
DO 16 K=1 .7 

16 UN I T ( I « J ) =UN 1 T ( I « J ) +SAVMAT ( I «K ) *PEMAT (K ♦ J > 

WRITE (6. 17) ( (UNIT( I* J> . J=1 .7) * 1 = 1 .7) 

17 FORMAT ( 9H I DENT 1TY/7 ( 7E16.7/ )// ) 
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15 WRITE (6.10) I KOUNT .PEVEC <7.1 ) * (PEVEC (I*1)»I=1«6) 

10 FORMAT (34H1 CORRECT IONS ON INITIAL COND I T I ONS ♦ 1 OX, 1 3H I TERAT I ON NO. « 
1 2X ♦ I 3///9H DELTA TF.4X.E16.8. 10X.10HDELTA P 

2SI 1 .2X.E16.8. 10X. 10HDELTA PS I 2 . 2X . E 1 6 . 8/39X • 1 OHDELTA PSI3.2X. 

3E 16*8. 1 OX ♦ 1 OHDELTA PS I 4 . 2X ♦ El 6. B/39X . 1 OHDELTA PS 1 5« 2X « El 6. 8 « 1 OX . 

4 1 OHDELT A PS 1 6. 2X« El 6.8// ) 

C 

TF»TF+PEVEC<7. 1 ) 

D051 I«1 *6 

SAVE < I +8 ) = SAVE ( I+8)+PEVEC( I . 1 ) 

51 CONTINUE 
C 

RETURN 

END 

BLOCK DATA 
COMMON /SPACE/ 


1 

VAR 

♦CUVAR 

♦ SAVE 

♦ C 

.MAXI T 

• IMAT 

2 

E 

• PE 

♦PEMAT 

♦PEVEC 

* ERRVAL 

♦SOMEG 

3 

DER 

• RS 

• BETA 

♦ DPI 

.TF 

♦ ELE1 

4 

RSV 

• FI 

♦ DP2 

♦ Cl 

« ELE2 

• DRSV 

5 

F2 

• El 

• I I 

♦SOMEGS 

.COMEGS 

• I PR I NT 

6 

F 3 

• E2 

♦ N 

♦COMEG 

♦ EN2 

♦KOUNT 

7 

Q 1 

* E3 

♦ TEMPO I 

• SIO 

»PH 10 

♦ Q2 

8 

E4 

♦ SPEC 

♦ STO 

♦LAMBDA 

.TEMPT 

♦ Q3 

9 

T2 

•TEMPSP 

♦ CIO 

• CRIT 

* ENl 

♦ KOUNT 

1 

RVXO 

• RVYO 

• RVZO 

♦DRVXO 

. DRVYO 

• DRVZO 

2 

REAL 

T 4 

LAMBDA ♦ MU 

• MU 

• CTO 

♦ B 

.XT 1 



DIMENSION 



1 

VAR (93 ) 

. CUVAR ( 93 ) 

. DER ( 93 ) 

2 

ELE1 <92 ) 

.ELE2 (92 ) 

.ERRVAL (92) 

3 

RSV < 3 ) 

.DRSV (3) 

. PX (6.6) 

4 

PPSI (6.6) 

. PPS 17(6) 

»CUPX (6.6 ) 

5 

CUPSI (6.6 > 

. CUPS 17(6) 

♦DRPX (6.6) 

6 

DRPSI (6.6) 

.DRPSI 7 (6) 

♦ ERX (6.6) 

7 

ERPSI (6.6 ) 

♦ERPS 17(6) 

* E ( 7 ) 

8 

PE (7.7 ) 

♦PEMAT (7.7) 

♦PEVEC (7. 1 ) 

9 

B (7 ) 

♦SAVE ( 93 ) 

* XT 1 ( 3 * 3 ) 


EQUIVALENCE 

1 (VAR (16) .PX <1*1)1 . (VAR (52 ) *PPS I < 1 . 1 > > • 
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2 (VAR ( 88 ) . PPS 17(1 ) ) 

3 ( CUVAR ( 52 ) . CUPS I ( 1 ♦ 1 ) ) 

4 (DER( 16) «DRPX( 1 * 1 > > 

5 (0ER(88) .DRPSI7( 1 ) ) 

6 (ERRVAUC51 ) .ERPSI ( 1 • 1 ) ) 

C 

DATA N/92/ 

DATA (SAVE < I ) « I =1 6.93 )/78*0.0/ 

DATA SAVE (52 )/l .O/.SAVE (59)/l . 0 / .SAVE ( 66 ) /\ * 0 / .SAVE ( 73 >/ 1 . 0 / * 

1 SAVE ( 80 )/ 1 • 0 / . SAVE ( 87 ) / 1 •0/ 

END 

SUBROUTINE I NT2 ( I I . N « NT * C I « SPEC « C I MAX « I ERR « VAR , CUVAR * DER »ELE1 ♦ 
lELE2.Et_T«ERRVAL«DERSUB»CHKSUB, I TEXT ) 

DIMENSION VAR (93 )♦ CUVAR (93 ) 

D I KENS I ON DER ( 93 ) * ELE 1 ( 92 ) « ELE2 (92 ) « ELT (92 ) « ERRVAL ( 92 ) 

DIMENSION TEMP (92 ) ♦ DER 1 (92 ) «DER2 (92 ) »DER3 (92 ) 

DIMENSION S 1 VAR (93 ) 

IF ( I I ) 1 * 1 *2 

C INITIALIZATION SECTION 

1 IF(CI) 3.4*3 

4 WRITE(6. 1000 ) 

1000 FORMAT ( 1 1 HOC 1 = 0 STOP) 

CALL EXIT 
C SAVE Cl 

3 H = C I 
18 I ERR= 1 

TO=SPEC+VAR( 1 ) 

MODE= 1 
11= 1 
N 1 =N+1 
DO 5 J= 1 *N1 
CUVAR(J)=VAR(J) 

5 CONTINUE 

C EVALUATION SECTION HERE 

8 CALL DERSUB 

IF (MODE.LE.l ) GO TO 6 

IF ( I 1-3)36.36.7 

36 CALL CHKSUB 
IF(II.EQ.2> GO TO 1 

37 DO 38 J=1 .N1 

38 V A R ( J ) = CUV AR ( J ) 

I F ( I 1-3)6. 7.7 

7 RETURN 

6 I F ( SPEC ) 9.7.9 

9 DEL=VAR ( 1 ) —TO 


* (CUVAR( 16) .CUPXt 1 . 1 ) ) ♦ 

, (CUVAR(88).CUPSI7(1 ) ) * 

. (DER(52 ) .DRPSI (l.l )) « 

* (ERRVAL (15). ERXd.l ) ) ♦ 


. (ERRVAL (87) «ERPSI7( l ) ) 
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DELP=DEL*(1 .+1 • OE— 6 ) 

IF ( ABS t DELP ) — ABS ( SPEC ) ) 2*10*10 

10 TO=VAR ( 1 ) 

GO TO 7 

2 11=1 

IF (MODE-4) 11.12*12 
RUNGE-KUTTA 

11 DO 20 J=2.N1 
DER3(J-1 )=DER2 ( J— 1 ) 

DER2 ( J— 1 >=DER1 (J-l ) 

DERI (J-l )=DER( J) 

ELE1 (J-l ) ssDER ( J ) 

CUVAR ( J ) =0 *0D+00 
DEl_T =0 • 4*ELE 1 (J-l )*H 
S1VAR(J)=VAR(J) 

CUVAR ( J ) =S 1 VAR ( J > +DELT 

20 CONTINUE 

S 1 VAR ( 1 )=VAR <1 ) 

CUVAR (1 ) =S 1 VAR ( 1 ) +0 *4*H 

CALL DERSUB 

IF( I 1-3 >23.23*7 

23 CUVAR ( 1 )=S1 VAR( 1 ) +0 .4557372S*H 
DO 24 J=2 *N1 

ELE2( J-l ) =DER C J ) 

DELT= (0.2969776 l*ELEl (J-l >+0. 1 5875964*ELE2 ( J- 1 ))*H 
CUVAR ( J ) =S 1 VAR ( J J+DELT 

24 CONTINUE 
CALL DERSUB 

IF ( I 1-3)25.25.7 

25 CUVAR ( 1 ) =S 1 VAR { 1 )+H 
DO 26 J=2 *N1 

TEMP (J-l )=DER( J) 

DELT=(0.21810040*ELE1 (J-l ) -3. 0509651 6*ELE2 ( J- 1 ) 
1+3. 832864 76* TEMP ( J-l ) >*H 
CUVAR ( J ) *S 1 VAR ( J ) +DELT 

26 CONTINUE 
CALL DERSUB 

IF( I 1-3)27.27,7 

27 DH=H 

CUVAR ( 1 ) =VAR ( 1 ) +DH 
DO 28 J=2.N1 

DOL'3=0. 1 7476028*ELE 1 (J-l )-0 .551 48066*ELE2 ( J- 1 > 

1+1 .20553560*TEMP( J-l ) +0.171 1 8478*DER ( J ) 

CUVAR ( J ) = VAR ( J ) +DH*DOUB 

28 CONTINUE 


oooonoo o on 


APPENDIX B 


MODE=MODE+l 
GO TO 8 
ADAMS— MOULTON 
ADAMS-BASHFORTH PREDICTOR 

12 CUVARU )=VAR(1 )+H 
DH=H/24.0 
DO 13 J*2*N1 

DOUB“55* 0*DER( J )— 59»0*DER 1 <J-l >+37«0*DER2< J-l )-9.0*DER3< J-l ) 

CUVAR< J)=VAR <J)+DH*DOUB 

13 CONTINUE 
DO 14 J = 1 » N 
DER3 ( J )=DER2(J> 

DER2 ( J )=DERl(J) 

14 DER 1 (J )=DER< J+I 1 
CALL DERSUB 
IF ( I 1-3 ) 15* 15«7 
ADAMS-MOULTON CORRECTOR 

15 DO 16 J = 2 « N 1 
TEMP=CUVAR ( J ) 

D0UB=9.0*DER(J)+19.*DER1 (J-l ) -5 • 0*DER2 C J-l )+DER3(J-l ) 

CUV AR ( J > =VAR ( J ) +DH*DOUB 

16 ERRVAL ( J— 1 )=(TEMP-CUVAR(J>)/14. 210526 

19 GO TO 8 
END 

SUBROUT INE MAT I NV ( A * N * B *M» DETERM « IPIV0T« INDEX«NMAX« I SCALE) 

*■**■»•*•***** REVISED 08/01/68 

MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS 

DIMENSION IPIVOT(N) « A < NMAX ♦ N ) « B ( NMAX *M ) ♦ INDEX (NMAX«2 ) 

EQUIVALENCE ( I ROW » JROW ) * ( I COLUM « JCOLUM ) « (AMAX* T« SWAP) 

INITIALIZATION 

5 ISCALE=0 

6 R1=I0. 0**100 

7 R2= 1 .0/R1 
10 DETERM=1.0 
15 DO 20 J=1 .N 

20 IP I VOT ( J ) =0 

30 DO 550 1=1 *N 

C SEARCH FOR PIVOT ELEMENT 

C 

40 AM AX=0 • 0 
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45 

50 

60 

70 

80 

85 

<50 

95 

100 

105 

106 


1 10 


DO 105 J=1 *N 

IF ( IPIVOT<J>-l > 60* 105* 60 

DO 100 <«1 »N 

IF ( IPI V0T(K )-l ) 80* 100* 740 

IF (ABS( AMAX )-ABS(A( J«K ) > )85« 100*100 

IROW=J 

ICOLUM=K 

AMAX=A( J*K> 

CONTINUE 

CONTINUE 

IF (AMAX) 110.106*110 

DETERM=0.0 

ISCALE=0 

GO TO 740 

IPI VOT ( I COLUM ) = I P I VOT ( IC0LUM)+1 


INTERCHANGE ROWS 


TO PUT PIVOT ELEMENT ON DIAGONAL 


130 IF < IROW-ICOLUM) 140* 260. 140 

140 DETERM =—DETERM 
150 DO 200 L=1 ♦ N 
160 SWAP=A ( IROW«L ) 

170 A ( I ROW « L ) = A ( I COLUM * L ) 

200 A ( ICOLUM.L >*SWAP 
205 IF ( M ) 260* 260. 210 
210 DO 250 L*1 « M 
220 SWAP»B< IROW.L) 

230 B ( I ROW « L ) =B { I COLUM * L ) 

250 B ( I COLUM « L > *SWAP 
260 I NDEX ( I « 1 ) * I ROW 
270 I NDEX ( I « 2 ) = 1 COLUM 
310 P I VOT* A ( I COLUM « I COLUM ) 

IF (PIVOT) 1000*106*1000 


SCALE THE DETERMINANT 


1000 

1005 

1010 


1020 


1030 


P I VOT I =PI VOT 

IF (ABS(OETERM)— R1 11030* 1010*1010 
DETERM =DETERM/R1 

i scale* iscale+i 

IF (ABS (DETERM) -Rl ) 1 060 « 1 020 . 1 020 
DETERM=DETERM/R 1 
I SCALE= ISCALE+1 
GO TO 1060 

IF ( ABS ( DETERM )-R2) 1040* 1040* 1060 
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1040 DETERM*DETERM*R1 
I SCALE* 1 SCALE- l 

I F ( ABS ( DETERM ) — R2 ) 1 050 « 1050.1060 
1050 DETERM*DETERM*R1 
I SCALE* l SCALE- 1 

1060 I F < ABS (PI VOT I l-Rl >1090* 1070*1070 
1070 P I VOT I *P 1 VOT I /R 1 
1SCALE*1SCALE+1 

IF <ABS< PIVOT! 1 — R1 ) 320. 1080. 1080 
1080 P I VOT 1 *P I VOT I /R 1 
ISCALE*1SCALE+1 
GO TO 320 

1090 I F ( ABS < P I VOT I ) -R2 ) 2000* 2000 « 320 
2000 P I VOT I *P I VOT I *R 1 
I SCALE* I SCALE-1 

I F ( ABS (PI VOT I ) — R2 1201 0« 2010 « 320 
2010 P I VOT I *P I VOT 1 *R 1 
I SCALE* I SCALE- 1 
320 DETERM*DETERM*PIVOTJ 
C 

C DIVIDE PIVOT ROW BY PIVOT ELEMENT 

C 

330 A < I COLUM « I COLUM ) * 1 • 0 
340 DO 350 L=1*N 

350 A(ICOLUM,L)*A< I COLUM. L > /PI VOT 
355 I F < M ) 380* 380. 360 
360 DO 370 L*1 »M 

370 B ( I COLUM «L 1 *B ( I COLUM, L ) /PI VOT 

C 

C REDUCE NON-PIVOT ROWS 

C 

380 DO 550 LI * 1 . N 

390 I F < L 1 — I COLUM ) 400. 550. 400 

400 T*A (LI . ICOLUM) 

420 A ( L 1 .1 COLUM ) *0 . 0 
430 DO 450 L* 1 «N 

450 A (L 1 «L ) *A (L 1 .L > —A ( I COLUM.L ) *T 
455 IF (M ) 550. 550. 460 
460 DO 500 L=1.M 

500 BILl «L )*8(L1 .D-B ( I COLUM.L >*T 
550 CONTINUE 
C 

C INTERCHANGE COLUMNS 

C 

600 DO 710 1*1 .N 
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610 L=N+ 1 — I 

620 IF ( 1NDEX(L» 1 >-INDEX(L«2> > 630« 710* 630 
630 JROW- INDEX (L«l ) 

640 JCOLUM= INDEX (L»2> 

650 DO 705 K=1«N 
660 SWAP=A (K« JROW) 

670 A(K« jROW)*A(K. JCOt-UM) 

700 A(Kt JCOLUM ) “SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 
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